{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.optimize import root"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.7473062499999997e-18\n"
     ]
    }
   ],
   "source": [
    "## constant\n",
    "hbar = 6.63e-34;\n",
    "me = 1.6e-31;\n",
    "\n",
    "print(hbar**2/(me*(1e-9)**2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEaCAYAAAAR0SDgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XuUHVWZNvDnobvp0LlgAk2CggkqChFC0BZsEbySKDqDzlIEE8FvdDUKo+CoBPBbg/O5BBpmBIRR6RlRNIiDzngZBwVkVFAaMFyGBIICSiSQhEBAAoGQdN7vj1NHTlftOteq2lV1nt9avdK1z+mqfXIu79l7v3tvmhlERER28l0BERHJBwUEEREBoIAgIiIBBQQREQGggCAiIgEFBBERAaCAICIiAQUEkZSQHCY5TvJXJK8k2ee7TiL1KCCIpGcNgLea2ZsA/AHA0Z7rI1KXAoJIm0j+D8n3x91uZo+Y2bPB4XYAO7KpmUh7FBCkVEj2kHyW5AGO264k+Y0Ozv1ikmuD35cDeAuAq0gaycPq/N0+AN4J4CftXjt0vmNJrib5DMkHSB6exHlFen1XQCRJZjZB8l4ArwawqlpOcgjAuwG8qoPTHwXgZ8HvnwZwDIABM9se9wckZwC4HMCHzOz5Dq5dPd+RAEYBfADArQD27PScIlVqIUghkPwoyWtIfpXkEyR/T3I+yVNI/onkYyT/Jrj7KgDzQ6c4H8A/mdkjda6xB8kfk9xA8imS/xV8oFcdBeDq4PeDANwH4HSS60g+SvKE0Pl6AVwJ4PNm9rsOHn6tfwTw/8zsZjPbYWYPm9nDwfX6g3pPkHw6+Jkg+RzJtyd0fSkxBQQpigUAhgB8H8DuAFYC+Glw28sBfAHA/w2O70alhQAAIPlXAPZFJSjUMwPAxQBeCmBecJ0Tg3P0ATgCwHXBfQ8K7rMawF4AzgZwVuh8xwE4FMA/kPwlyQ+EL0jyJySfjPn5Sei+PcH/wSDJ+0muJXkJyV0AwMy2AjgFwPVmNs3MpgF4AMAiM/t5g8cuoi4jKYyDAJxjZtcDAMl7APSb2UXB8Sq88HpeBeD4oLwHwLkAPmdmW4KyL6Ly4b4BwPHVcjO7H8D9wTm2krwOwMzg+AgA/2tmm4PjBQC+bmb/EZxzBSoDx39hZt8G8O16D8rM3t3C/8FsAH0A3gfgcADbAPwIlUD4ueA+B6ISLEFyKirBclXkTCIOaiFIUSzA5EHZ+Y7je4Pf7wawb/Ct/iMAngPwLQAIBptfbmaHA/g5gL+tnoDk+0n+Juj+eRLA6QB+H9xc210EVALUz2qO9625flqqGUsXm9k6M3sMwJeCulX9JSAEv68zs00p10tKQgFBco/kXAA744UPZwBYCODOmuMFNccPohIEDgbweQB/by/sBHU4Xuhq+imANwbXeCsqg7WnAngxKt1Fj9ac8ygA/x3ctw/AfqHrh+vT7GP7aU1/f/jnp7X3NbMnAKwFUG9Xq9qAcFDN7yINqctIiuAgACvNbAfwl8yduQDuCt3nhwBgZkZyNYCvAbjFzH5Vc7+ZANYFv/8ZwKyav38IlW/5MwH8M4A9ANwTpI32m1m1BTAfwJOhAeqFAL7c6gMzs3e2+CffAPAJkj9DpcvoVAQtJZK7AxgEcE9w3/0wOYiK1KUWghTBQYi2Bu6vGRPYCcABofusCspOC53rCQC7Br/vCqDanXIFKv3z61H5gL0PwD1Bqui7MLm7aAGA/w2dt7aFkqYvAPgtKh/0qwHcAeCLNXV4oGYy3MMAjiF5aAb1khKg9lSWbkLyQABnmNkHSY6g8s3/4gZ/czWAS8zs6nr3Eyk6dRlJVzGzlSTXkLwRlTGC45v4s18C+EWqFRPJAbUQREQEgMYQREQkoIAgIiIAFBBERCRQqEHl3Xff3ebNm+e7GiIihXLbbbc9ZmaDje5XqIAwb948rFixwnc1REQKheSaZu6nLiMREQGggCAiIgEFBBERAaCAICIiAQUEEREBoIAgIiIBBQSRDixbBpAv/Eyd6rtGIu0r1DwEkTyZMQPYvHly2ZYtlcCwZAmwfLmfeom0Sy0EkTZMmRINBrWuuAJYvDi7+ogkQQFBpEWLFwNbtza+37XXAuPj6ddHJCkKCCItuvba5u97xBHp1UMkaQoIIi2IW1txp5h30vbtlYFnkSJIPSCQvIzkoyRX1ZTNInkdyfuCf2emXQ+RTo2PA2scS4T19QETE0Dc5oPnn59uvUSSkkUL4ZsA3hEqOx3A9Wa2L4Drg2ORXDvmGHf588+/8PuSJdHbzdRKkGJIPSCY2Q0ANoWKjwZwefD75QDek3Y9RDq1dm207JBDJh8vX15JOw370pfSqZNIknyNIcw2s3UAEPy7R9wdSY6QXEFyxcaNGzOroEituBTSW26Jln32s9Gy7duVcST5l/tBZTMbM7MhMxsaHGy44Y9IKq67LloWbh1UjY66WwknnJBsnUSS5isgbCC5JwAE/z7qqR4iDY2PuweMXa2Dqg9+MFp2333J1UkkDb4Cwo8BVL8vnQDgR57qIdLQ8cdHy6ZMqf83cctWjI11Xh+RtGSRdnolgHEAryK5luRHAJwL4EiS9wE4MjgWyaUHHoiWffKTjf9u2rRo2VlndV4fkbSkvridmR0Xc9Pb0r62SKfiuotGRxv/7UknAeedN7nsUXWOSo7lflBZxKePfzxa1uwS166gsWOHso0kvxQQROq4++5o2cknN//3rm6jk05qvz4iaVJAEKlj+/ZoWTPdRVWuD/+VK9uvj0iaFBBEYriWm+jra+0cruAxMdFefUTSpoAgEuMrX4mWxa1nVM8uu0TLtLaR5JECgkiMLVuiZe1si/m610XLLr209fOIpE0BQcRhfLySEVSrt80k7XMds2yeeaa9c4mkSQFBxOF0x4Lszaabhg0PRzfQcQ1Wi/imgCDicPvt0bITT2z/fAMD0TKNI0jeKCCIONRuegNUvuG3km4a5ko/1TiC5I0CgohDODW03fGDqtHRaLfRs892dk6RpCkgiISMjUUDws47d37ecFDRfATJGwUEkZBzzomWveY1nZ83HFQmJrQctuSLAoJIyIYN0TJX6mirXEHl7LM7P69IUhQQREK2bZt83NtbSR3tlCuoaJtwyRMFBJEa4+PROQJJjB8AlaASHkfQfATJEwUEkRquCWmDg8mdXxPUJM8UEERquCaknXlmcucPT1DbsUMT1CQ/FBBEaoTXL+rrA0ZGkju/61yXXZbc+UU6oYAgUiM8N2DKlGTPPzoa3VNh69ZkryHSLgUEkcDYWPTDuacn+euEA4ImqEleeA0IJD9F8m6Sq0heSTLh72MizbvoomjZggXJX8ds8rFaCJIX3gICyZcA+CSAITM7AEAPgGN91Udk/fpoWRIT0sJmz558PDGhgWXJB99dRr0AdiHZC2AAwCOe6yNdLLzY3M47JzMhLeyMM6JlWvlU8sBbQDCzhwH8E4A/AVgH4M9mdm34fiRHSK4guWKjpnVKilwzlNMwMhI9d3i5bREffHYZzQRwNIB9ALwYwFSSS8P3M7MxMxsys6HBJGcIidQYG4tOEpsxI73r7bLL5OO0go9IK3x2Gb0dwB/NbKOZbQPwnwDe4LE+0sVcA8qvf31213/uueyuJRLHZ0D4E4DXkxwgSQBvA7DaY32ki7kGlE87Lb3r7bbb5ONt2zSwLP75HEO4BcD3AdwOYGVQF60OL16Exw8GBtIZUK5yDSxrxrL45jXLyMzOMrP9zOwAM/uQmSkjW7wITw5LY0JarZERoL9/cpnmI4hvvtNORbwbHwe2bJlclnZAAJJbVlskKQoI0vVcS16nMUO5EaWeim8KCNL1Vq6MlqUxQzksnHq6dav2WBa/FBCk65GTj6dPT3dAuerDH46WXXhh+tcViaOAIBISXo00LaOj0YHlDRuyubaIiwKCdL2nnpp8HE5BTVN4YDnLa4uEKSBIV1u2LLpkxaxZfuoCaG8E8UsBQbraN78ZLUtyD+VGwmsYbdlSSYMV8UEBQbpauHUwMJDsHsqNHHhgtOy887K7vkgtBQSRGknvodyIK7315puzrYNIlQKCdLXwpjhZD+oOD0fnIzzzTLZ1EKlSQJCuNT4eDQhTp2Zfj3BAyCrtVSRMAUG6lquvPss9EKqUeip5oYAgXcvVV5/mHghxwuMWmzdrCQvxQwFBula4rz7tPRDiLFwYLdMSFuKDAoJ0rXBffdYZRlWuVskTT2RfDxEFBOla4b57X/sTDA/7nR0tUqWAIF0rvIaRT9osR/JAAUG60rJl0V3SZs70UxcgujlOnoKVdA8FBOlKrjWMTj0182r8hWtNI2UaSdYUEKQr+V7DKEyb5UgeKCBIVwr32c+Y4aceVaOjlaBUS5lGkjWvAYHki0h+n+S9JFeT9JAFLt0ojxva+w5KIr2N75KqiwD8zMzeR3JnAAON/kCkU+PjwKZNvmsRFQ5SeQxaUm7eAgLJGQCOAPBhADCz5wHoLSCpc61h5DPDqCq8htFzz/mph3Qvn11GLwOwEcA3SN5B8t9IRtaaJDlCcgXJFRs3bsy+llI6d94ZLfOZYVQVXmlVu6dJ1nwGhF4ArwHwVTM7GMAzAE4P38nMxsxsyMyGBgcHs66jlFD4m/f06X4zjKpcK61q9zTJks+AsBbAWjO7JTj+PioBQiRTPvZAcHGtaXTHHdnXQ7qXt4BgZusBPETyVUHR2wDc46s+Ir651jTautVPXaQ7+c4y+gSAK4IMoz8A+D+e6yNdIM9bVGpNI/HJa0AwszsBDPmsg3SXsbHKBjS18pBhFEepp5IlzVSWrnLOOdGyPGQYxdm0SZlGkh0FBOkqjz8++bi/Px8ZRlWufRGUaSRZUUCQrhLeJS0vGUZVp5wSLVOmkWRFAUG6Sl52SYszMlKZF1FLmUaSFQUE6Sp5zjCqylurRbqHAoJ0jaJlGIlkTQFBusZFF0XL8pxhVFWEVo2UgwKCdI3wktd5WcMoLJxptHmzttOUbCggSNfKa1+9K9NI22lKFhQQRHLGlWmk7TQlCwoI0jWK1Bef19aLlJsCgnQFV4ZRf7+fujRD22mKDwoI0hVcGUYHH5x9PZoV3k6zSK0bKS4FBOkK4QwjwL0hTV7sttvk461blWkk6VNAkK40a1ZlQ5q8OuOMaJkyjSRtCgjSlfK2hlGYMo3EBwUE6QpFHJRVppFkreWAQHImSaZRGZE0jI9HxxDynGEk4kvdgEDyH0juF/zeT/IXAB4AsIHk27OooEinXBvM5DnDKE4RWzlSLI1aCB8A8Lvg9xOCfwcBvAnA2WlVSiRJd94ZLctzhlHVlCmTj7WdpqStUUB43sws+H0xgO+a2YSZrQbQm27VRJLx3HOTj/OeYVS1cGG0TNtpSpoaBYStJA8gOQjgLQCurbltIL1qiaQn7xlGVa5WjLbTlDQ1CginAvg+gHsBXGBmfwQAkkcBSOSlSbKH5B0kf5LE+UTKYng4uhS2ttOUNNXt9jGzmwHs5yi/GsDVCdXhFACrAcxI6HwikxR52YeitGakHOoGBJJ/X+92M/tSJxcnuReAdwH4IoC61xJpR9m2zVSmkaSpUZfR9Jqfz4SOp9f5u2ZdCOA0ADvi7kByhOQKkis2btyYwCWlm5xzTrSsCNtmxlGmkaSJLyQRNbgjeYeZJZa9TfLdAI4ys5NIvhnAZ8zs3fX+ZmhoyFasWJFUFaQLzJgxuYXQ3x/NOsqzV78auOeeyWXveQ/wgx/4qY8UE8nbzGyo0f1amancXORo3mEA/prkgwC+C+CtJJcnfA3pcn19k4+LthyEaztNZRpJWrytZWRmZ5jZXmY2D8CxAP7HzJb6qo+UU3hQtmiDtK5F7pRpJGlpNKi8Ei+0DF5B8q7qTQDMzBakWTmRTpVhEHbq1OjAuEgaGs02rtunnxQz+yWAX2ZxLekeWtROpDWN5iGscZWT7EGlm8d5u0gelGVRu7Aiz6uQfGu02ukMkmeQvITkIlZ8AsAfAByTTRVF2lPURe3CwrOVN2/WdpqSjkaDyt8G8CoAKwF8FJW1jN4H4GgzOzrluol0pKiL2oW5Mo20naakodEYwsvM7EAAIPlvAB4D8FIz0xCXFE7RMoyqRkaAz3xm8sCyttOUNDRqIWyr/mJmEwD+qGAgkr2izZ+QYmrUQjiI5FPB7wSwS3BcTTvVgnSSWxp8FWlNoyyjnqwqIpKksi1qF1aG+RWSP95mKoukqWyL2oVpkTtJgwKClNLjj08+7u+vDM4WVTj1FNB2mpI8BQQppaIvahemRe4kCwoIUkpFX9QuTIvcSRYUEKSUyphhVPRWjuSfAoKUTtkzjETSooAgpXPRRdGyMmUYVZWxFSR+KSBI6YSXvJ4+vdgZRlVa5E7SpoAgpROetBXOOCoqV6bR2WdnXw8pLwUEKZ1t2+ofF9XISDRbKtwaEumEAoKUzsTE5OOytBAAYNq0ycdlemzinwKClMrYGLBly+Sy2bP91CUNZZtfIfmigCCl0i0ZRlXKNJIkKSBIqZQ1w6hKmUaSJm8BgeTeJH9BcjXJu0k6cihEWlPWDKMqZRpJmny2ELYD+LSZ7Q/g9QBOJjnfY32kBMqaYVSlTCNJk7eAYGbrzOz24PfNAFYDeImv+kg5lDnDqEqZRpKWXIwhkJwH4GAAt/itiRRZ2TOMqpRpJGnxHhBITgPwHwBONbOnHLePkFxBcsXGjRuzr6AURrdlGFU9FXnXiLTHa0Ag2YdKMLjCzP7TdR8zGzOzITMbGhwczLaCUijr108+HhgoV4ZRVTjTaMsWYNkyP3WRcvGZZUQAXwew2sy+5KseUh7hAeSeHj/1SJsr0+iyy7Kvh5SPzxbCYQA+BOCtJO8Mfo7yWJ/SmjIFIKM/O+1Urhz28OBqWQdbR0YqrZ9aO3b4qUsaDj3U/XolgT339F27cvOZZfRrM6OZLTCzhcHP1b7qU0bz51feRHFbLZoBJ54IzJiRbb3S8txzk4/LPNgafs7K8lhJ4NZb429fv75yn6VLs6tTN/E+qCzp6OkBVq9u7r6bN1feZEW2bFk0w6jMu6SFJ+AVfQmL8fHWXoNXXKHWQhoUEEqIbK8LochB4ZvfjJaVOcPIbPLx5s2VD9UiGh8H3vCG1v9u/Xpgt92Sr083U0AomU4HUova9bB9++TjsmYYVR14YLTsvPOyr0cS2gkGVZs2AfPmJVaVrqeAUCIzZtRvGZi98LNTzDO/bVsx32DhQFaWcZE4554bLbvjjuzr0al6X2DmzHnh9bpkSfz91qxR2m1SFBBKYvHiSreBy6xZ0S6GiQlg7lz3/desKV72UbgPPdzHXjbDw5WVXGvFPf95NW9e/BeYm24C1q174Xj58uhruFZRW0d5o4BQAuPjwLXXum+bOxd4/HH3bQ8+GB8UTjwxkaplYnw8+mFYpjTMZhVpIb+xscoXD5ebbqoEPJd6rduyzjvJkgJCCbzxje7yWbMqH/r1PPhgNKe96tBDO6lVdlzfDhcsyL4eWQvPswgv7JdnH/+4u3zJkvhgUBX3OHfsqKRaS/sUEApu8WL3t+GenviWQVhcymK9fPA8ufnmaJmrj71swuMkW7YUo6tv6VL3a3bWrErXUDPiuo9Wry5utlUedEVAWLy48gHZ11e+CS1xXUXhrJtGLr3UXV6EXO9wQBsYaPwtswwWLoyWXXhh9vVo1RVXuMub/QJTddpp7vIjjmjtPHk3b14lJXzmzPSDHa3eSE3OzJq7vx15ZmuLttx1F/DEE5PLBgaA170uwYp5Mj7uHjzde2/gZS9L7nwHH5zvrJ1f/3pyN0JPT3w3Wpk89VQ0s6i3FzjsMD/1aYbr/Qi0/5r9zW/cX372mA3sv1/r58ubG26ItobaeT9e9bE33GZmQ43uV/oWwp//HC3bsgVYfW/2dUnSunXuD2+yvTcWEP+tetWq9s6XlW4cQAYqHwrhAda8f79zBYNOXrNxwe/RDe2dL09uvNH9fD75ZIoXNbPC/Lz2ta+1Vi1aVJt9P/nnpptaPl1u9Pam85ji/r/y6tJLo3WdO9d3rbIza9bkxz59uu8axVuyJJ3XbNx5Fy1Kpt4+zJ2b7OcWgBXWxGds6VsI11wTv+plJzMkfVq2zN1M3n//zvvOr7nGXZ7X7I1zzomWnXlm9vXwJTwhL89LWLjGDvr7O3/NLl/uTkWNG1/Lu2XL4lNyDzkk3fGx0gcEoP4kpSIu1XD++e7ye+5J5vyHHBIta3ahvKyFByL7+8u9ZEVYeLMcIJ+TtOKSOb785WTO/9WvusuLkjpdK+75mz4duCXlTYa7IiAA8X2r27ZVspCKYulS92OpN7W/VXEvuiK8uYoY4Dvh2izHlYbr25VXRsv6+pIL3iMjwLRp0fKipE5XTZniLu/ry2ar1K4JCEBlBqRLkZqWcSl7zeZvN8vVSsjjm6vV9NqyGRmJBsG8LYU9NuYe+L/kkmSvE/c+LsraXEuXxu9dktVSLF0VEIaHgUWL3LcV4ZtlXEsmLh+7E3GthDzN4xgfB559dnJZWXdJq6e/33cN6vvsZ6NlPT3Jd+0ND7uXYlmzJr/jKrXivuxlmTnWVQEBqAyaut5AReg6cn0D6ukBRkfTuZ6rf/o730nnWu1w9bXOnp19PfImbwv7ubo6Pv3pdK4Vt1TLUTnfnDeuqyjuC2xaui4gANGtFqvy3HUUF6y+8pX0runK4DHLz/IIrr7yMm+KE2fq1MnHW7fm5zmKG3dK60sM4O7ufPLJ/LYS4rqK+vris/7SUqiZykNDQ7ZixYpEzrV0qbuJNn16NoM3rXLtZjZtWvpLHvf0RPt/Z81qfZmBNAwMTO4y6u0t1oqfSXnve4Ef/nBy2dy5jRc2zILrdbtoUfofdL7eL+2I26kwyY9mkpqpXM/y5e5VPjdvzt9mG3HfsrJo0Rx3XLRs06b0r9uM8KqXRRgHSoNrDCkPz1Hc+yiLb72urpann85Py6kqbq2wJLMGW9G1LYSqLKJzp1x1fNGL3MsAZHX9JUuSz2xqhWsf3ry0XHyYOrWyJEtVHlq6U6ZEu0KyfI5cr9v+/vgu46yNjbn3HUmjpVuIFgLJd5D8Hcn7SZ7uow5xgzZ5mZkb1zq4+urs6vCKV0TL4jIisnK649WiAeUX+P7QGx9394u7xqXS4npv52l85WMfc5ffcEO29ajlLSCQ7AHwLwDeCWA+gONIZv4xfM017m8SeVlX3ZX7PziY7fLO3/qWu9znG2vlymhZNw4oV4W7y7Zt8/v8HH+8uzzLWeRxXVOf+lR2dYizeLG7FyKJ5Wc64bOFcAiA+83sD2b2PIDvAjjaR0W+9jV3+ZvfnGk1IuIyi370o2zrMTzszu8/44xs61ErPCGt25asCHM99rPPzr4eVfffHy3LOoUy7ppbtvj/shc3/pfU8jPt8hkQXgLgoZrjtUHZJCRHSK4guWLjxo2pVGRkBJgzJ1r+/PN+B5hdLxpfm7+4vlX5HLgMd0fE7bPbLUZHo3sKp/R2acjnYHKz1/Q5LyFu5nQaE0xb5fNt5BrOjTSizGzMzIbMbGhwcDC1yqxb5y73tVBYXOvggguyrUdVXN64j5nLro18tMF6tNvIV2LExRdHy1yTHLOSp3kJ4+PulUynT093bkazfAaEtQD2rjneC8AjnuoCID7Vy8cMZlfrYNo0v90irsHlq67Kvh6uAeXXvCb7euRNb+/kY1/rPIWXEwGyHUwOi1uGxUcXVlw3tO+MsCqfAeG3APYluQ/JnQEcC+DHHuuD5cvdA8xZz2D2Oe+gHtfgso+JYLffHi0799zs65E3u+wy+Xjbtuy7POOu53t8Jw/zEpYtcy8rsv/+2dWhEa/zEEgeBeBCAD0ALjOzL9a7fxrzEMKWLXN3E82ZE9+tlDTf8w7q2WmnaFdE1nMS+vsnv7G6dYZymOu1u+uuKW+5GBKePQ5k+96px/W+6unJriXlc85TIeYhmNnVZvZKM3t5o2CQldHR6DctAFi/PptvE3EDTlnOO6jn5S+PlrnWuk9T+A3c7QPKVaOj0f8LV/dNmlzX+8d/zLYOcVythImJbFpRcfOafM1IjqO3ksP117vLXbMKkxQ34DR3rt/c5FqubqMdO7IboFu2LLq2kiuAd6vwOEJ4eY80xY21+e4uqorLOEo7cWR83L3jIOl3tr+LAoLD8HB8v16au4bFDTjlYZGyquFhd0bPCSdkc31XKy3tQF0k4UyjiYns+smvuy5a5mPgth4fiSOHH+4uj5v/5JMCQoy4CSK33prOt+G4ASdXypxvxx4bLbvvvmyuXbteD1DpIslDul5euLKtspigNjbm7gv3MfegnqwTRxYvdrfS5szJT8uplgJCHXETReIififimq1pb6rdjrhmbhbfRMNvrnAXSbdzZVtlMUHtrLOiZXGbvvjm2sENiF95tBNxgSYPg+wuCgh1jI5WJoyETUwk23VUlAGnWq6JRmkvZTE2piWvGxkejgbJuH16k7R+fbTsk59M/7rtyCpxJO61mef3tQJCA3ETRpLqOirSgFMt10SjtJeycH0L1YS0qPC6U2mPI8Rl6eS5Ky/txJFDD3WnQg8M5Pt9rYDQhLiIHl6Pvx2HHeYuz+OAU624/s80l7J47LFomSakRbmWAU9zHCFvS1U0o17iSFzqd7PGx92rFAPAM890du60KSA0IW53NSC6n20r9tzTPRA3a1Y+B5zCXEtZpDknIdxd1NOTn3TcPHF13aU5jpC3pSqaFZc4smZNZy2quC+Kee4qqlJAaFJcZN+ypb3xhKVL3f2uQHF2/cpyTsKyZdHgqfkHbiMj0QlqaW2Yk/e5B43EJY6023UUN27Q35/vrqKqQm2hOWvu/nbkmZd5u/66dcDvf+++bY/ZwP77NXeep54C7rjDfdveewMve1l79fPhhhvcH9RJp8v+5jfRGcq77gosXJjsdcrC9by88pXJZ9L86lfRspkzgQULkr1Omlyvrao3van589x4Y3TSZDvnScNVH3tD/peuKJo994xPpXt0A3DXXc2dJy4Y9PYWKxgAwOAe0bI0lktw5XIX7f8qSwOOrsykJzjGpU4WKRgA8eN4gDvgudQLBq98Zet18qVQLYQsFrdrRm8ICzwlAAAI+klEQVRv/JIA9TYRj1s4D6hkFcW9oPLONdHntNOSzTIJX2OnnbJdlqFoxsejfdlJLwK4227RzLIpU7JfPykJrv+vWvU+JuMWrQMqA9e+d0EDCrK4XVHVWx1x06bKCyScbdPbW3/NlKIGA8CdUXLhhcmd35W5pPGD+oaHo+MISa/q6Uozzuvcg0aGh+sP+pLR7rb58+sHg4GBfASDViggtKlRw+qKKyovlupPvW+zBWqkObkySp5/PrnB5e99L1p28snJnLvMXAOcSaUFxw0m53nuQSPLl9ffm2D9+snvadf8oaq+vvynmLqoy6hD9b4hNKNA//11ufZJ2Hff+EH4Vrj+j8vy/5ampUsrX0xq9fW518xqVZrPt2/z5rlXHW5Wf396WV3tUpdRRszaX4+/TB9qRx4ZLUtiwTtXPrj2T26OK80xiTGE8XH3a/fyyzs/dx48+GD7q7TOnZu/YNAKBYQETEy0lmY5MFCuYADEr2rZaReFa7mKwcHOztlNXMGz02UsjjkmWrbTTuWaJHjNNcBNN7X2N5demq+l6tuhgJCQW26pfMjXm7Lf01N5kRWxb7EZrsfe6czlDRuiZXnZgasIXMHTFWRbsXZttOy44zo7Zx4ND1fe041aC4ccUrlfUSbj1aOAkLDHH6+8OFw/27eX61tUmGtwuZOZy3Fr7JfhjZcVV/CMmyHfjLjB5CLMwm3XNdfEv6fN8rlEfbs0qCyJcg027rUX8NBDrZ/Llec+dSrw9NPt168bJTlPxHWusgwml5kGlcUL1+Cyq4uhGa48d6Wbts7VlfflL7d+nrjxoLIMJosCgiQsbnA5bhOgOHEfPkXOc/fF1ZXXTiZMOIUVqCRIlLkbtNt4CQgkzyd5L8m7SP6A5It81EPS4cq4qjeJx8U1GL3vvu3Vp9uNjLi7elrZWD4uQF9wQXt1knzy1UK4DsABZrYAwO8BpLz5omQpbpCt2ZU2x8bcS3moa6J9rq68VjaWd7UO+vo0wF82XgKCmV1rZtWVVW4GsJePekh6XK2EZvesPeWUaFnZ8tyz1sk8kbj9Pi65pP36SD55zzIi+V8A/t3MnIlrJEcAjADAS1/60teu6WROuWTK1U3RaFXXsTH35iRLlpQ7tTELu+7q3iO80UeA63lMagkMyYb3LCOSPye5yvFzdM19PgdgOwBHg7TCzMbMbMjMhgY1RbVQXK0Es/o7zJ10krtcwaBz55/vLq83lhC3A5haB+XkrYVA8gQAHwPwNjPb0szfaB5C8cQt/nfTTdEuoMWL3f3aixbFd3lIa/r73d/sXR8D8+e7kwHmzInfHEfyyXsLoR6S7wCwDMBfNxsMpJji1pgPb0YyPh4/yKlgkJyLL3aXhwP30qXxmWEKBuXlK8voEgDTAVxH8k6SX/NUD0nZ8uWVXHWX6kZCS5fG71ZVb9MSad3ISPx6W9WgsNtu7qwiQM9H2XkfVG6FuoyKq519I/K4rnxZ6PnoLrnuMpLu0+pSwoA+fNJ02mmt/42ej/JTQJBMDA+39iHUTgCR5o2O1t8uMqxAHQnSAQUEyczoaGUTkUZcGUiSvHvuaTwm4Fq9VspLAUEyNTJS+YCZMyd625w5ldsUDLKzfHnl/7y/P3rbokWV3QCle/T6roB0J6Uu5ovGBwRQC0FERAIKCCIiAkABQUREAgoIIiICQAFBREQCCggiIgKgYGsZkdwIoIg75OwO4DHflchYtz3mbnu8QPc95iI/3rlm1nBDmUIFhKIiuaKZhaXKpNsec7c9XqD7HnM3PF51GYmICAAFBBERCSggZGPMdwU86LbH3G2PF+i+x1z6x6sxBBERAaAWgoiIBBQQREQEgAJCqki+n+TdJHeQHArddgbJ+0n+juRiX3VME8nPk3yY5J3Bz1G+65QGku8Insf7SZ7uuz5ZIPkgyZXB81q6jc5JXkbyUZKraspmkbyO5H3BvzN91jENCgjpWgXgbwDcUFtIcj6AYwG8GsA7AHyFZE/21cvEBWa2MPi52ndlkhY8b/8C4J0A5gM4Lnh+u8Fbgue1jLn530TlvVnrdADXm9m+AK4PjktFASFFZrbazH7nuOloAN81s61m9kcA9wM4JNvaSUIOAXC/mf3BzJ4H8F1Unl8pMDO7AcCmUPHRAC4Pfr8cwHsyrVQGFBD8eAmAh2qO1wZlZfR3JO8KmuCla2Kju57LWgbgWpK3kRzxXZmMzDazdQAQ/LuH5/okTltodojkzwE4dgjG58zsR3F/5igrZP5vvccP4KsAvoDKY/sCgH8G8LfZ1S4TpXkuW3SYmT1Ccg8A15G8N/hWLQWmgNAhM3t7G3+2FsDeNcd7AXgkmRplq9nHT/JfAfwk5er4UJrnshVm9kjw76Mkf4BK11nZA8IGknua2TqSewJ41HeFkqYuIz9+DOBYkv0k9wGwL4BbPdcpccGbpuq9qAyyl81vAexLch+SO6OSLPBjz3VKFcmpJKdXfwewCOV8bsN+DOCE4PcTAMT1ABSWWggpIvleABcDGATw3yTvNLPFZnY3yasA3ANgO4CTzWzCZ11Tch7Jhah0oTwI4ES/1UmemW0n+XcArgHQA+AyM7vbc7XSNhvAD0gClc+Q75jZz/xWKVkkrwTwZgC7k1wL4CwA5wK4iuRHAPwJwPv91TAdWrpCREQAqMtIREQCCggiIgJAAUFERAIKCCIiAkABQUREAgoIIi0i+XTN70cFq1++1GedRJKgeQgibSL5NlTmmSwysz/5ro9Ip9RCEGkDycMB/CuAd5nZA0HZX5G8heQdJH9OcrbfWoq0RhPTRFpEchuAzQDebGZ31ZTPBPCkmRnJjwLY38w+7aueIq1Sl5FI67YBuAnARwCcUlO+F4B/D9Zw2hnAHz3UTaRt6jISad0OAMcAeB3JM2vKLwZwiZkdiMq6TVN8VE6kXWohiLTBzLaQfDeAG0luMLOvA9gVwMPBXU6I/2uRfFJAEGmTmW0i+Q4AN5B8DMDnAXyP5MMAbgawj8/6ibRKg8oiIgJAYwgiIhJQQBAREQAKCCIiElBAEBERAAoIIiISUEAQEREACggiIhL4/3dAeu00NwvSAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1bd2d04bdd8>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.3306690738754696e-16\n"
     ]
    }
   ],
   "source": [
    "P = 2*(3*np.pi/2); # larger P gives us more solution...\n",
    "KA_scan = np.linspace(-4*np.pi, 4*np.pi, 4000);\n",
    "\n",
    "plt.figure();\n",
    "plt.plot(KA_scan, np.cos(KA_scan)+(P/(KA_scan))*np.sin(KA_scan), '.b')\n",
    "plt.axhline(1);\n",
    "plt.axhline(-1);\n",
    "plt.xlabel('Ka');\n",
    "plt.ylabel('RHS')\n",
    "plt.title('$mV_0a/\\hbar^2 = 6\\pi$')\n",
    "plt.savefig('solving_kp_transcendental.png', dpi = 300)\n",
    "plt.show();\n",
    "\n",
    "def RHS(x):\n",
    "    return np.cos(x)+(P/(x))*np.sin(x);\n",
    "\n",
    "## roots at pi\n",
    "print(RHS(np.pi)+1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Notes on the Transcendental Eq.\n",
    "The solutions are values for which the  value of the root func is less than 1. Cuz then we can solve the left hand side."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [],
   "source": [
    "def RHS(x):\n",
    "    return np.cos(x)+(P/(x))*np.sin(x);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAU8AAAFFCAYAAABsVm+UAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xl8nHW59/HPlaRpkjZ0oS1LC1SEFpG1RIqyngdUPILIARGQc1j04ehx4zwcHwQUBBW3o+L+iCzCUREQDygeUVCRRbYUurC3YAuFSlMKtGn25Hr+uHI7k3TSZKbJzD2T7/v1ymsm99yT/JrOfOe33+buiIhIfqpKXQARkXKk8BQRKYDCU0SkAApPEZECKDxFRAqg8BQRKYDCU0SkAApPEZECKDxFRApQU+oCFGrGjBk+d+7cUhdDRCrMokWL1rn7zOHOK9vwnDt3Ls3NzaUuhohUGDNbNZLz1GwXESmAwlNEpAAKTxGRAig8RUQKUNTwNLOrzWytmT026PjHzexpM3vczL5azDKJiBSi2DXPHwNHZx8ws38AjgP2cfc3A/9Z5DKJiOStqOHp7ncD6wcd/gjwZXfv7D9nbTHLJCJSiDT0ec4DDjWzB83sz2b2lqFONLOzzazZzJpbWlqKWEQRkYHSEJ41wDTgIOBTwI1mZrlOdPcr3L3J3Ztmzhx2AYCIyJhJwwqj1cAvPa5E95CZ9QEzAFUtpWDusHgx9PXB8uXwvvfB0qXx2H77Qe6PZ5GRS0N43gL8L+AuM5sH1ALrSlskKQfu0NwMX/0qzJgBO+4YoTh3LrzwAnz3u9DeDhs2wKpV8J3vxHM+9jHYaSdYuRJefDGOVVXBWWfBAQcoWGVkihqeZnY9cAQww8xWAxcDVwNX909f6gJOd10PWYbQ3Q2nngo9PRF8jz4a9wczg8sug7e/PVPzfMc74He/gwsuiMAc7Ec/gsMOi9rqHntE2NakoXohqWTlmlNNTU2ujUHGh+5uOPnkqCm+8AIMHivcc88IveyaZ3U1vP/9UaPM1tcHN9wAvb0Da5733ANPPDHw3JkzYc4cOPDAqMUqSMcHM1vk7k3DnqfwlDRyhwcfhNNPhzVrYOPGgY8vWABTp8KRR8J550VYbo2eHjjnnPhdTzwBTz018PFp02DWrOhHravbut8l6TbS8NRnqaRKTw989KNw++3w/PMDH9t5Z9h9dzjzTDjllM1rlVujpiZqlxC105/9LJrtK1fC2rXw6qvxNWVKhOgHPwif/ezWh7aUL9U8JRV6euAjH4Gf/jQGebK9+c1Ru/zAB0Y3MEeitxe+9CW4/vrNm/VTpsBJJ8H3v68mfSVRs13KQnd3DOb8/vcDQ7O2FubPh4ceSk8zuaMj+j+XLRt4vK4OTjsNfvADhWglGGl4pmGSvIxD3d3w3vdCYyPcemsmOKur4cQTobU15mWmJTghyrJ0KXR2xih+Ugvu6IArr4x/y//+37lH/6XyKDylqHp74eKLYfLkCM3OzjiehGZ7O9x0E0yYUNpybkltbaamfNRRmePZIXrNNdF3KpVL4SlF4Q733QeTJsGll0JXVxw3K5/QHKy2Fu64Iz4ABofoWWfFB8SmTaUrn4wthaeMufZ22G47OOSQTE0TInA6OsovNAdLQrSjAxYuzBxvb48A3W23eEwqi8JTxkxXFxx6KDQ0DJzYvtdeESx33BHBUykmToQHHoC2NnjjGzPHn302/gbveU/09UplUHjKqOvrg6uugvp6uPfezPHGxhgIWrYsXQNBo62+HlasiIn9yb/THX7968zfpEwnuUgWhaeMqra2mP/4oQ9lBkyqqmIA5bXXos9zvJg8OT4sLrkkc6y3N2rj2223+XxWKS8KTxkVSRN90qQIjMSJJ0Z/3xlnFH+CexpUV8NFF20+qNTSEk35s87S1KZyNQ5fzjKa3KMZOriJPm1ajDSX+2DQaEkGlVpbIzQT11wTTfvsDxwpDwpPKVhHR6w3P/TQgXMar7kG1q0bGBISJk2K/UUHN+UbG2GffQbORpB0U3hK3np74TOfidrm6tWZ4wsXju8m+kglTfn29ph5kEgG0q68UhPsy4HWtkte2tpi1/bB69BfeSUGSCR/GzdGN0dvb+bYtGmx12h9fenKNV5pbbuMqq6umOQ+adLA4Pzwh6NvU8FZuMbG+JueeGLm2KuvakAp7VTzlGFt3AjbbDPw2LRp0WRXv+bo2rQpavbZK5ImTIgwHU/TvEpJNU/Zap2dMYgxODg1IDR2Jk2KD6sPfzhzrLs7avbHHqsVSmmimqfktGFDTHbPtnAh/PnPsQxRxl57O+yyy8ClrVVV8Prr6iYZS6p5SkE6OmJddnZwVlVFmD7wgIKzmOrr4W9/i5p+oq9P05rSQuEpQLwpr7wy3rDPPZc5fuSRUQNqbCxd2cazqqqY+tXeDrvumjmeTGvasKFkRRv3dNEAoaMjmodr12aOVVfHIIVCMx3q6mKzkb/8JRYlJL1tU6bA3nvDww+rVVBsRa15mtnVZrbWzB7L8dh/mJmb2Yxilmk8y57snh2cyebECs50MYODD44PuyOPzBxPaqF3363dmoqp2DXPHwPfBa7LPmhmOwFvB57P8RwZA8kGxdnXQ6+pidqmBiPSrbYW7rwz/u+mTs2sRjr8cJgzB5Yvr+wt/9KiqDVPd78bWJ/joW8C/xfQ5+YY6+2FCy+MaUbZwXniibF6SMFZPpLJ9dm10NWroyXxox9piedYK3mfp5m9B3jR3ZeYWamLU9FyLa2sq4s5m5qAXZ6SWmhrayxcSFYjnX12XOteSzzHTklH282sAbgQuGiE559tZs1m1tySPflNtqinB848c/Ollf/6r1H7VHCWv8mT48PxhBMyx5Ilnp/5zMB18zI6ij5J3szmAre5+15mtjfwB6Ct/+E5wEvAge7+ty39HE2SH5nW1ugXy37zNDbG/EGtEKpMmzbB9OmZK5RC1D61KmxkymKSvLsvc/dZ7j7X3ecCq4EFwwWnDK+7O5bzNTYODM7PfS5TI5HKlOzmn73RSHt7HNdGI6OnqDVPM7seOAKYAbwMXOzuV2U9vhJocvd1w/0s1TyHtnFjzP/L/q+dORNWrVL/13izaRPMmhVN+kRNzfi7nlQ+UlnzdPdT3H0Hd5/g7nOyg7P/8bkjCU7JLXsjj+zgvOaaaKYrOMefXDvX9/REH+kxx2ijka2hjUEqQHIdocMOG3h8113h8cc1509CW1vMA3311cwxs9hoRAsiMlJZ85TRl0x2Hxyc994by/kUnJJoaIhBo6uvzhxzj5bK3ntro5F8KTzLVLK0sqFh4JZlyXWEDj44ahUi2aqqYtra4OsnPfaYrp+ULzXby1CuqShaWimF2LAhprJlx8B4v36Smu0VKJl+NHnywODU0kop1DbbREvlqKMyx3T9pJFRzbNM5Jp+pMnuMppaW2HbbQd+MFdXx7Sm8fTBrJpnhUh2dh88/eiSSzTZXUbX5MnRJZR9/aTe3viQPvjggaEqCs/Uyt5rM3tn9732is7+iy6KWoHIaKqpgR/8ILqBZs/OHP/LX2Kz5Xvu0Z6hiZLvqiSbyzUgpJ3dpZjq6+GFF+C++2Ln+sRhh2m1WkI1zxRJVgjlGhDSzu5SbGZwyCHRdbRwYeZ4S0t0Fx1yyPhuyis8U6C3Fz772fgkX7Ysc3zatKiF3nQTTJhQuvLJ+DZxYlw5dePGgYsu7rtvfDflFZ4l1toaL8gvfGHz9ejaQkzSZPLkeL1mr5OHaMong03jicKzRDo6YLfdoimePZfuqKOi+X7GGbEaRCRNqqtjsHJwUz6ZZ7zvvuNnmafenkXW0wMf/GA00Z99NnN85sz45L7jjri0gkiaJU351taBffFLl8Zr+7Ofrfzd6xWeRdLXB1ddFU307I0Zqquj7+jll9VEl/IzaVJMor/mmsxeCu7RDVVbO/Aig5VG4VkESb/mhz408NM4GUV/29u0iYeUr6qq6GYavMyzry8Wd2y77cDNmCuFwnMMtbfHROPGxoGbziYT3TWKLpWktja6ndraYlVcYv36qKHutlsEbKVQeI6Bjg7Yffd4wbz0Uub4nDnRr7lsmfbZlMpVXx97yW7cOLA/9Nln47FDD62M+aEKz1HU2Qn77RehuWJFZupRXV28kF54Qf2aMn5Mnhy71N9338AW1r33xvvgsMPKO0QVnqOgqyteCA0NsGRJZjPZCRNiTlxr6/jalUYkYRZ9+m1tm284cs898Z457rjyvJaSwnMrdHXB4YdHMN5zTyY0q6rihdLWpg08RCCz4UhX18BLIvf2wq9+Fe+h9763vEJU4VmAzk7Yf//4D7/77sx/uFm8MDo64oVSo21XRAaYMCEGSjs7B47Md3XBrbfGe+rww8ujOa/wzENHB8ybF//BixdnQrO6OkKzs1Mj6CIjkYzM5wrRu++O99j++6d7tZLCcxjusZfhrFkxcrh8eWY55YQJmbmaCk2R/A0O0WRJcnd3VFAaG+O9l8Z5okUNTzO72szWmtljWce+ZmZPmdlSM/tvM5tazDINpbsbjj8ett8+dtFuacmEZmNjpk9ToSmy9ZIQ7eiICkkyla+7O957jY0x1e/ss9NzXaVi1zx/DBw96NgdwF7uvg/wDHB+kcv0d319cO21MUdzyhS45RZYuzbz+BveEI+/9pr6NEXGQtIn2toKl14Kc+fG8b6+uKLnj34U+0DsvDPcf39pt8Ir+gXgzGwucJu775XjseOBE939A8P9nNG6AJw7PPQQfOQjcTG1NWsGPl5bC/Pnxzma2C5SfO3tsUn4c89tfk35+fNjoPbRR0fv/VmuF4A7C/jtUA+a2dlm1mxmzS0tLQX/kr4+uO46WLAA9twTDjoo/vhJcFZVxSfeRz4SK4KWLlVwipRKfX2MNXR0wEknwR57ZB57+ml46qnYOHzBAjjgAHjwweLUSFNT8zSzC4Em4J98BIXKp+bpDs3N8OUvwyuvRB/KE08MPGfixAjMnXeG227TtnAiadbZCUccEbfZC1MSBxwAO+0UM2HOOw+amka++c5Ia56p6LUzs9OBY4AjRxKcQ3GHRx6B22+PybfJuvJ166L/cvD+glOnwq67xq4vCkyR8jFxYvR5wsAgXb48+ksXLYoviPf+4YfHxiRmsOOOMV5x3nlbt4Cl5OFpZkcD5wGHu/uIJyS0tcWnTfKp8/TTsHIlfPGLQ09r2GWXqF1WVUVz/fLLNegjUu6yg7S7G047LQZ6p06N6U4rV8If/xhf2ZYuhWOPjUCdPz9yYd99ASaN6LqgRW22m9n1wBHADOBl4GJidH0i8Er/aQ+4+4dz/oAsEyc2+XXXNXPuuRGWr72W6ed43/ti27fsHY1qahSWIuNNTw+cc06MZ8yYkal5Pv443HhjnGMWQdvQAF//Opx88j6d7kuHHeUoep/naNlzzyZ/7LHmATVP9/j0eP/7df0fERlaXx/ccEPcDq55VldPfsK99c3D/YyyDc/RmqokIpKtXKcqiYiUBYWniEgBFJ4iIgVQeIqIFEDhKSJSAIWniEgBFJ4iIgVQeIqIFEDhKSJSAIWniEgBFJ4iIgVQeIqIFEDhKSJSAIWniEgBFJ4iIgVQeIqIFEDhKSJSAIWniEgBFJ4iIgVQeIqIFEDhKSJSAIWniEgBFJ4iIgUoania2dVmttbMHss6Nt3M7jCz5f2304pZJhGRQhS75vlj4OhBxz4N/MHddwf+0P+9iEiqFTU83f1uYP2gw8cB1/bfvxZ4bzHLJCJSiDT0eW7n7msA+m9nDXWimZ1tZs1m1tzS0lK0AoqIDJaG8Bwxd7/C3ZvcvWnmzJmlLo6IjGNpCM+XzWwHgP7btSUuj4jIsNIQnr8CTu+/fzpwawnLIiIyInmFp5lN2ZpfZmbXA/cD881stZl9EPgy8HYzWw68vf97EZFUq8nz/JfM7OfA/3P3h/P9Ze5+yhAPHZnvzxIRKaV8m+1fI2qHD5jZo/2j35PHoFwiIqmWV3i6++eAucDxwEvA94na6A/MbL9RL52ISErlPWDk7n3u/it3fzfwRuBbwHuARWb2oJmdYWYTR7ugIiJpsrWj7RuIFUOtgAFTgKuAFWZ2yFb+bBGR1CooPM3sYDO7DngRuAT4I7Cvu+8BvAl4DvjhqJVSRCRl8hptN7OPA/9KBOSTwKeA69x9Y3KOuz9jZhcTm3yIiFSkfKcq/SdwC/BRd//zFs5bDlxacKlERFIu3/Dc2d1fHu4kd0+a8yIiFSnfqUrDBqeIyHiQb5/nH7fwcB/wOrAIuEpBKyKVLN9muwHzgB2AvwIvA9sBbwDW9H//j8C/m9nh7v7EKJZVRCQ18g3PbwCXAwe4+6PJQTM7ALiR6OdcBPwe+CKxEmlMtLVBXx8sWRK3y5fDSSdBVRr2iRKR1HOHxYvjNmE28ufnG55fAD6XHZxRCF9kZpcAX3D3vc3sa8TI/Jh59lm46SY499wI0g0b4L77YOFCWLky/iAvvQQ77ABLl8KsWfDtb0NNvv9iESlb3d3wgQ/A9OkRjDvumHns5ZfhxhsHnl9bCzCpfiQ/O98omQesG+KxFmC3/vvPApPy/Nl5eeMb4X3vg3nzouZ57bXwne/Ad7879HPuuAPmzIGpU+G112CPPeI5ClSR8tfdDaeeCj098f6eNg2am+GFF4Z+TmMj/PCHMH9+fG8GCxZsah/J78s3NlYCHwJ+m+Oxs/sfB5gBvJLnz85LQ0M00fffP77ff39461uhtzd3zfORR2DFivhK3HUX3HwzvOlN8Sl0yinwmc9AdfVYllxERkNPD3zsY/DQQxF6L7wAuS5ttuOOcOyxm9c8q6rg6KNhwYL8musJ8+wG/3Anm50C/AR4HLiZuGTGLOAEYC/gVHf/uZl9H9je3f8p/yKNTFNTkzc3N4/4/O5uOO00WLs2ap5PPRVfg82eDZMmxR9z8WKoqxvFQotIwdzhwQfhjDOgvh5WrYJXX938vAULYJttouY5YQL85CdxO1Jmtsjdm4Y7L6+ap7tfb2briIGhC4AJQDfQDLzD3e/sP/X/AL35/OyxNmEC3HBD5vu+PvjZz6LZXlcXNdPWVnjxxcw5U6ZEX6kZPP10/IeJSPG4wwMPwPvfH+/D55/f/JxZs6I7bto0OPPMaEEWY+A4r5rngCeaVRHN83Xu3jeqpRqBfGuew+npgXPOiU+2DRvgmWcGPl5TA9tvryAVKYa2Npg7N95va3NcEnK//aCrK4Ly/PNHt6tt1GueZlYL/A04o38/zz4q6EqXNTUDB5s6OuDAA6NZsHp1hOvq1fFYY2ME6aRJMVVKTXuRrdfRAfvsE++59eujdZht993jvfbQQ+l4z404PN29y8x6gI4xLE9q1NXFQBNkgnTNGli3Lgalkub95MnRh3r88fCDH2jkXiQf3d0xa+bee2OEvHdQZ9+cOfF1110wMWVbrOfbM3ALcOJYFCTNkiBduxbuvz+mSTU2xmO9vfDKK3DllXFsyhTYtKm05RVJs74+uOYamDEjKh+33hrvoSQ4t9km3mNtbTGCfv/96QtOyH+q0m+Bb5vZL4ggXQMM6DR19y2tfy9rZnDQQTHdqa8P/uu/4IILYkoURA21oyNeENtuC29+c8wtjYm3IuNbZ2fm/dPauvnjs2fHwpeDDips6lCx5TtVaaiBISfWvbu7F2WW5GgPGG2Nri5497sztdNs1dUxJ3XNmugjFRlP3OEvf4F3vStaZIP7MbfbDg49NGa+5DOdaCyNyVQl4B8KLM+wzOzfiQn4DiwDznT3suhfra2NGiZEU2P+/MzgUm8vbNwYtdHGxhglvPNO1UalsnV2xlLpFSs278Yyg113hWXLynvWSr7zPLe0e3zBzGw28AlgT3dvN7MbgZOBH4/F7xtLDQ3RT5MsFfvDHzITeTduhHvuiRfM1KlxXkNDacsrMppaW2HmzGiNDa5lTp8em/dUypLoQi8AN8PMjjGz081sev+xuv65n4WqAerNrAZoIK4LX7YmTIj+m/Xr45N3zpzMY319cXzSpKiNfuhDMRVKpBx1d8Nxx8VrubEx+v2zgzMZ/HnllcqakZJX2Fn4GrAa+BVwNTC3/+FbgQsLKUT/ZTv+E3ieGIR63d1/X8jPSqOkNtrVBSeeOHCOWmsrXHVVHJs9G9pHtCWBSOm1tcXAaF0d/OpXAweBpk+HD384gnXFivJung8l35ri+cDHiIu7LSQGiRK/Bo4ppBBmNg04jthUeUdgkpmdluO8s82s2cyaW3LtAJBySW20rS060WfOzDzW2xuj9g0N0T969dWbN3tESq23Fz772QjMSZM2n8w+a1aEaKXVMnPJNzw/BFzq7pcBjwx6bAXwxgLLcRTwV3dvcfdu4JfA2waf5O5XuHuTuzfNzE6eMmMWO0CtXRtNnIULBy4v27QJPvjBeHEedljUWEVKqbMT9t03QvMLX4jvExMnwlFHxbGXXx4/s0ryDc/ZwANDPNZF4Xt4Pg8cZGYNZmbAkcR14SvexImx8UFXV0wcnjIl81hHRwwwNTTExGFNvpdico8Nxhsbo9m9dOnAvvmZM6MF1d4+Pucz5xueLxJbz+WyL3Fdo7y5+4PAL4ja7LL+cl1RyM8qV1VVsdXWa69tPsCUPd2pri5qpRpgkrHS3Q3veU98sB9ySDTDs6eD77VXBObatdGCKocJ7WMh3/C8CbjIzA7OOuZmNg84F/h5oQVx94vdfQ9338vd/9ndO4d/VmVKBpg6O6M5lN0M6uyM/tC6Othtt6idioyG9vYYtKyvh1//OkI00diYGQBatiwdG3OUWr7h+TngKeBuYHn/sZuI2uJy4MujVjL5++T7jRtzDzA9+2y80DXAJIXq7YWLLorXUENDDFpmb86RDABt2FD5A0D5yns/TzOrBk4F3knsIv8KcDvwU3cvWmMyTcszi6mzEw4/PLblGvxfV1cXq5sefDCdGylIenR0wN57x+bCgwcka2rgiCPgN78Zf/2YMPLlmXlPanf3Xnf/L3c/zd3f4e6nuPu1xQzO8SwZYOruhksuiYGkREdH7C/a0BA10lybL8j45R5bv9XXx2tkxYqBwTllSgxadnaOzwGgfOkq52WqujqaW6+/nnsFU0dH9FPV1cXGC5ruNH51dMRGwpMnx2uho2NgqyVZAfTaazFoWYxLWFSCfFcY1ZrZxWb2lJm1mVnvoC/VPktg8AqmadMyj3V2Rm0jme50332bN/el8vT1RT94Y2OmltnWlnl8PKwAGmv5dv9+Dfgosa/nL4FxOyKeRskKJoja6K67ZrbIS6Y7HXJIvFGamrS7UyXq7IyrHjz55MDR8sScOXENLm1Is/XyDc8TgYvd/YtjURgZPZMmxWqPXLs7tbdnJt/X18O3vqXmWjnr7Y3+729+M2qXlb6bUVrkuxny68DxadgtfryOtm+NwbXRbHV1scnD8uVqwpWLTZtgxx3jAzLXhjJvfGP575lZCmM12v5r4LDCiiSlltRGc/WNdnTERe2SbfK0pj6dOjtjQ+3GxhgA2rBhYHCqL7N48q3Efwe4rv9yHP8DrB98grs/NxoFk7GT3Tfa3h7z/V58MTMK29qaadZPmQL/9E+aIF1KyRUm77knZlcMvsLkjBlxKeyHH9bKn2LammsY5XzieLyGUSVwj/mjxx0HuXb7q6uL0L38cvWPFkPSj3n55VHbzNUKmDULnntu/OxiVCwjbbbnG55nMERoJtz92hH/wK2g8Bw7PT3w8Y/DjTfGfo2DTZsW80y/+lU4/XQF6Wjp7YVLL4Xvfz8Cc+PGzc/Zdlv4xCfgwgsHbmMoo2dMwnOYX1gNTHH3HG+30afwLI4tXRkUFKRbaySBmcYrTFayUQtPM1sPHOXuj/R/b8QlN87J7t80s4XAX9Rsr1wdHTGHcP366CMdbPr0GKB4y1ui1qo3em5dXfCOd2SuX/7665ufM2dOfN11l/YpKLbRvPTw1EHnVRGX2/hcYUWTclVXFzVQyB2kSRP/xRejeTllSgw6LVmigYz2dthjj1iUsHp17q0EFZjlReOnUpDsIO3sjF14XnwR/va3GB3euDHTBG1shLlzY2rNtGlw++2Vv7Kpqwve+c7MB8tf/7r5Btb19bHN4Pz5cNttlf83qTQKT9lqEyfC/ffH/Z4e+OQn4U9/ivBcvTqOrViROX/aNJg3L3Ygnzq1MsK0sxMOPjhW97jDM88MXEue2G236Oc880y44AIN+pQzhaeMqpoa+N73Mt+3t8M++8T9jo4I07Y2WLw4c87228Mb3hA1UzN405vSvZSwuxtOOQXWrYvv3WHRotzXmJozJ2qYdXWxB+t4776oJCN9ec42s13771dnHXst65w5iAxSXx9LPiFqZT/9KXzlKzGY1NYWNbRXX82suwf4859jP8kkdBMLFsD55xevtpbUont6MnNf3WNX/1wzDyZMiOv7uMdzrrwy+oXH6zV+Kt1IRtv72Hxupw11TKPtMlLu0Nwczdf166Pm+fLLsevPUA46KC6Bu+OOEUpmsMsusSP63LnRf7h8edzut1/0y+67bzx3yZII5MWL43fsvnvcrlwJO+8Mq1ZFwK9ZAzvsAL/4RawNH8qCBdGf6x59nBroqQyjOdp+5iiUR2QzZjGt6Y47Msf6+uD66+GWWwYuQ3z11aiRPvBAfA2lsTGaz9tuG03/Cy6Am2+Ox044AS67LBYAvPJKrMzJNa9ysOOPz9xPlq9edlls66da5fg1apPki001z/El6Ve85poI1WLUPF9+OQZ2DjhAITmeFH2FUbEpPEVkLIzZBeBERCRF4WlmU83sF/3XR3rSzN5a6jKJiAwlTTPpvgXc7u4nmlktoKusiEhqpSI8zWwbYof6MwDcvQvQPuYiklppabbvCrQA15jZo2Z2pZlttsWrmZ1tZs1m1tySa8deEZEiSUt41gALgB+4+/7AJuDTg09y9yvcvcndm2bOnFnsMoqI/F1awnM1sNrdH+z//hdEmIqIpFIqwtPd/wa8YGbz+w8dCTxRwiKJiGxRKgaM+n0c+Gn/SPtzaFmoiKRYasLT3RcDw87qFxFJg1Q020VEyo3CU0SkAApPEZECKDxFRAqg8BQRKYDCU0SkAApPEZFh6c2TAAARpklEQVQCKDxFRAqg8BQRKYDCU0SkAApPEZECKDxFRAqg8BQRKYDCU0SkAApPEZECKDxFRAqg8BQRKYDCU0SkAApPEZECKDxFRAqg8BQRKYDCU0SkAApPEZECKDxFRAqQqvA0s2oze9TMbit1WUREtiRV4Ql8Eniy1IUQERlOasLTzOYA7wauLHVZRESGU1PqAmS5HPi/QONQJ5jZ2cDZADvvvHORiiVp4g6PPAK/+x3stBM8/zzssguYxde8ebBiBbzvfbBsGeyzDyxZEs/dd9+47575WU8/HbfusHJl3FZXw3nnxa3IUFIRnmZ2DLDW3ReZ2RFDnefuVwBXADQ1NXmRiidF0tsLX/oSPPoozJgB69bBzJkDz1m/Hn7zG2hr2/z5ZjB5cjy2ahVccQVcdhmce248/vWvx/2urvi+pwdeey0TptnuugsOPRReeinzeEtLplwAxx8Pp54KValpv0kxmed65RS7EGZfAv4Z6AHqgG2AX7r7aUM9p6mpyZubm4tUQhktfX3wk5/At74FjY0ReFOnRoi1tMDjj4/s55x0EhxzzNjUPO+5B+64Y2TlePObYffdB4awGZx1FnzgAwrWcmRmi9y9adjz0hCe2fprnv/h7sds6TyFZ7q5w0MPwb/9W9xPaoWrV8Nzz235uW96U9T6ctU8q6rg4IPhlFPGLph6e+ErX4nboWqeTzwBTz01/L+jvj7z/KoquO8+mDhxbMoto2Ok4ZmKZruUt54e+NjH4MEHo8nc0ACvvhpN56FMmxa1wuya57bbwoIFpe9vrK6GCy7Y8jlJwP7hD7DNNgNrnkuXxr//yRzzRrbdNmrHmzbF36mtLZr+n/mM+ljLTepqniOlmmdpdHfD+98Pjz2WObZuXYRFLnV1UQPL7o9cuBAuvxxqKvSju7sbTjsN1q6FDRsyofrYY/FYLjNmxAdK4h/+Ab73vcr9G6VZ2TbbR0rhOfZ6euCjH4U//SkCoKMDXn8dNm7Mff6UKbDddlGjam+HnXeG226D2triljutOjvhiCPiNql5PvlkfJ/LtGkRqu3t8SF02mmqoRaDwlPy0tsLn/88XHVVvFE7OqC1NZqjuUyZArNmxf2amqgpfetbqinlK6mlLl0a/wcAL76YezYBRJjW1WX+jz7/efiXf9HA1GhSeMqQ3KN/8l/+Jd6kZvFGTKbgDFZbG7XIjo7o3zvlFDj/fNWAxkpPD3zyk1Hj7+mJmudLL8VMhVy23TbCtL4+pmE99VTcl8JowEj+rqMD9tsvajbt7VHbWbt26PN32y2eM3lyDOr85CcwYULxyjve1dREf2e2rq7M1KxNmyIsV6yIx155ZeC5jY2w/fYRoB0dMWPh/vs1yj/aFJ4VpqsL3vGOeGOZxZts5cqoweQye3acN2mS+ijTrLYWfv/7gceSQH366fhgrK+P//fe3mj6J1avjhZD0npwh7e8BW68UR+KW0PhWca6u2OyeHNzpg9s/fqh+8tmz45bsxiMeOiheJ6Up1yB2t4erYX29gjTlpYY5OvqytRUIcJ16tR4HSQ11He9C77/ffVbj5T6PMtEVxe8853xoh/JyHd9faYvzD2a3gsXRnDK+JGs6LrkkhhU6uiIZn57e+7zJ0+OJn9HR7x29t4bbrhhfNVQNWBUppLBnNNPH7jmeksjsBrQkXz09saa/2uuiQ/l+vroS03W/A/W0AA77pj54K2qgsWLK7fVovAsA93dcPLJsVxx06Y41tkZL+ShTJkSAwAdHTEwoClCMhp6euATn4Df/jZeSx0dMcF/w4bc59fWxp4CEP3lbW0xh7USJvYrPFOkqwuOPjozZ9I9XmwtLUOvzDGLUe/EhAmqUUpxZddQq6szNc/ly4d+zsyZ0Y/a0JA5/w1vgJ//vHya/grPIkteaHfeGZPLJ0/O9EeuWDF03yTEC27q1LhvBtOnx5ZomloiadTRAQceGLcQNc/lyzOtp1ymTYO5c+N+Y2O8Rxob4YMfTN/uUwrPUdbXBz/7Gfz3f8f9V1/NbGqxfj08++zA6SG5zJsXL7Sk5jllChx3HHz606pNSnnr6YFzzon++g0bMjXPdevghRe2/Nw99xy4e9a0adFKmz493htNTcUd6FR4jlBvL3z5y5v/B7tHs9oslsQ9+ig8/PDwP2+nneKFkF3zrK2Fb3875tZptFvGk6TS8Z3vZDZFSWqeq1dvebEGRKXi8MMzff3r1mW2BZwxI6bfrVkDc+aMXiVkXITnww83s2RJZsPbXJvcJtxji7Sdd47bZAPc3/4WHnhg5L93l11g//03r3m2t6sWKZKPpLb6xBMDl54mNc+//nXL2xoO9va3xz6wCbNMV8Hzz8fjSX/t/PkDKzJmsQovNtUeB+F55ZXNnHBC5lILI7m8Qi4HHhiBmG1wzXPdOthhh8reSk0kTZJwXbMm3o9D1TyXLoVbbx16FV0i2e0LotKT/T6urY3VdRGg42Bt+777ws03R81z3rzCap5r1yoQRdKopga++93hz0suCnj77QNrsPnWPPfdN7/ylXXNM02j7SJSGUZa80zRBAERkfKh8BQRKYDCU0SkAApPEZECKDxFRAqg8BQRKYDCU0SkAKkITzPbycz+ZGZPmtnjZvbJUpdJRGRL0rKupgc4190fMbNGYJGZ3eHuT5S6YCIiuaSi5unua9z9kf77G4EngdmlLZWIyNBSEZ7ZzGwusD/wYI7HzjazZjNrbmlpKXbRRET+LlXhaWaTgZuBc9x9s6unuPsV7t7k7k0zs3dPFREpstSEp5lNIILzp+7+y1KXR0RkS1IRnmZmwFXAk+7+jVKXR0RkOKkIT+Bg4J+B/2Vmi/u//rHUhRIRGUoqpiq5+72Aru4jImUjLTVPEZGyovAUESmAwlNEpAAKTxGRAig8RUQKoPAUESmAwlNEpAAKTxGRAig8RUQKkIoVRoVavBj23TfuL1kC++wTt+5xzB2eeQbmz4f99wfTGiaRcccdHn0Unn4a5s3L5IBZ5MfSpZEdS5dm8mQkyjY829rghBPg5pvj+xNOgMsug3PPha6uONbTA6+/DttuC//+77DLLvHHe+YZ6OuDVaviGMT9nXeO2yR8s730EuywA6xZk7l1hx13zNwH2Gkn+PSnobp67P8GIuWqtxe+/GV44YXMMbN4b730Uub+mjXxHsvFLN6/zz8f793nn4/vs8Nx3jy48074xjfglVdgyhSo6U+92lr4+tfhggsiOy64IJMnI1G24dnQAD/7WeaT4uab49Nj3rzNa54rV8KFF8axyZOhtTV3QI6WG26At74185/Y0gIzZsC6dUPfJuWZOTO+339/OP98hbCkS28vfOlL8Mgj8b1ZvIZbWjL3c72+Z86M13hyf8kSeOCBsS2rWbzfN22CL3wB5s7NXfOcPz+yY/78/Gqe5mOZImOoqanJm5ubR3RuX18EmvvY1jxbWuCXo7gT6V57wdveNjBcIROwgwN38PcAxx8Pp54KVerdHteSmt7zz8frFAYGX67vcx1bsQKWLRudMlVVwXvfG6/X5HeNRc1zxQo46aSRvwfMbJG7Nw173ngIz2Jxh+Zm+NrXYNq0wmuejzwSP2e0LFwYn6jJ75k5M8qWlGvw/S2pqoKzzoIDDlAf8ljp64Prr4f77oPtt48wySW7JpfrfnaN8KWXRremd8QRmdd4ITXP6mo488x0vo4UnmUsefPccgtMn154zXPVKli0aPTLV1cHxx4bNePsmvgOOwz9RkhqCatWxfdz58ZtUvsf/Lyk1rB8Oey+e9zOmxePPfNM3K+qig+FZJDQDPbbL85ZsiQeG/xz3TODi8lAQfbz99kHbropfqdZputn993jFjKtl+zuoVytmKRGNFRrxj1Ts3LP1LjWr4cbbxzxf0dektCDwmqeNTXwqU9BU1P6Qm+0KDxlyBDemprnihXwhz+Mfdmz+6smTYrbKVPisddfj/sNDdHhnwwS1tbCbbfFOclgYhKmicWLM4OLyUBB9vM//vHoH99mmwiKZNBx0qToK4fi9JsDHHkkHHbY6NQ8t9sODj4YTjlFXTjDUXjKmBjcpFTNc/RrnrNnw667wsknK+hKQeEpIlKAkYanPtdERAqg8BQRKYDCU0SkAApPEZECKDxFRAqQmvA0s6PN7GkzW2Fmny51eUREtiQV4Wlm1cD3gHcBewKnmNmepS2ViMjQUhGewIHACnd/zt27gJ8Dx5W4TCIiQ0pLeM4Gsnb2Y3X/MRGRVEpLeOZa1LfZ0iczO9vMms2suSV79wIRkSJLy2bIq4Gdsr6fA2y2HYK7XwFcAWBmLWa2KsfPmgGsG4tCjqFyK7PKO/bKrczlVl4Yusy7jOTJqVjbbmY1wDPAkcCLwMPAqe7+eAE/q3kk61LTpNzKrPKOvXIrc7mVF7a+zKmoebp7j5l9DPgdUA1cXUhwiogUSyrCE8Dd/wf4n1KXQ0RkJNIyYDSarih1AQpQbmVWecdeuZW53MoLW1nmVPR5ioiUm0qseYqIjLmKDE8z+7yZLTWzxWb2ezMb4uKl6WBmXzOzp/rL/N9mNrXUZRqOmb3PzB43sz4zS+0oa7ntmWBmV5vZWjN7rNRlGQkz28nM/mRmT/a/Hj5Z6jJtiZnVmdlDZrakv7yXFPyzKrHZbmbbuPuG/vufAPZ09w+XuFhDMrN3AH/sn3XwFQB3P6/ExdoiM3sT0Af8EPgPd0/dNVH690x4Bng7MZf4YeAUd3+ipAXbAjM7DGgFrnP3vUpdnuGY2Q7ADu7+iJk1AouA96b1b2xmBkxy91YzmwDcC3zS3fO+MHNF1jyT4Ow3iRyrldLE3X/v7j393z5ALBJINXd/0t2fLnU5hlF2eya4+93A+lKXY6TcfY27P9J/fyPwJCleWu2h/zqoTOj/KigfKjI8Aczsi2b2AvAB4KJSlycPZwG/LXUhKoT2TCgiM5sL7A88WNqSbJmZVZvZYmAtcIe7F1Tesg1PM7vTzB7L8XUcgLtf6O47AT8FPlba0g5f3v5zLgR6iDKX3EjKnHIj2jNBtp6ZTQZuBs4Z1PJLHXfvdff9iBbegWZWUPdIaibJ58vdjxrhqT8DfgNcPIbFGdZw5TWz04FjgCM9JR3RefyN02pEeybI1unvO7wZ+Km7/7LU5Rkpd3/NzO4CjgbyHqAr25rnlpjZ7lnfvgd4qlRlGQkzOxo4D3iPu7eVujwV5GFgdzN7g5nVAicDvypxmSpK/wDMVcCT7v6NUpdnOGY2M5nNYmb1wFEUmA+VOtp+MzCfGA1eBXzY3V8sbamGZmYrgInAK/2HHkjz7AAAMzse+A4wE3gNWOzu7yxtqTZnZv8IXE5mz4QvlrhIW2Rm1wNHEDv+vAxc7O5XlbRQW2BmhwD3AMuI9xvABf3LrVPHzPYBriVeD1XAje5+aUE/qxLDU0RkrFVks11EZKwpPEVECqDwFBEpgMJTRKQACk8RkQIoPEVECqDwFBEpgMJTRKQACk+pCGZ2cv+G0p39m9web2Z39a9dTjbB/Wb/xiatZvY3M/u1me1R4qJLmSrbjUFEEmZ2FJkNYM4llox+i9irMdlzdCLQCHwBWANMB/4NeMDM9nD3vxW73FLetDxTyp6Z3QdMA/Zy977+YwuJjaX/7O5H5HhONRGoLwMXufs3i1diqQRqtktZ6w/BtwC/SIIToH+D25WDzj3JzB40s9eIfVM3AZOJTWRE8qLwlHI3g2iev5zjsb8fM7NjgRuIy0ScCiwkQrcFqBv7YkqlUZ+nlLt1QDewXY7HtiO2JITYy3OFu5+RPNi/ie/0sS6gVCbVPKWsuXsvsenxiWb299dzf5/n3KxTG4imerZ/JvZ1FMmbap5SCS4Gfg/cYmY/JEbbLwGyR9BvB95rZt8EbgMOAD5BbOQskjfVPKXsufudxFVS5wO/BD4FnENmmhLAj4AvAu8Hfg28GzgWeL2ohZWKoalKUrGSCfK5piqJbC3VPEVECqDwFBEpgJrtIiIFUM1TRKQACk8RkQIoPEVECqDwFBEpgMJTRKQACk8RkQL8fxPFdZcNlmLDAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1bd2d7235c0>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "## do a scan of K...\n",
    "Kguesses = np.linspace(1e-3,4*np.pi, 10000);\n",
    "band_structure = [];\n",
    "for kguess in Kguesses:\n",
    "    val = RHS(kguess);\n",
    "    if(abs(val) <1):\n",
    "        q = np.arccos(val);\n",
    "        E = kguess**2;\n",
    "        band_structure.append([q,E]);\n",
    "band_structure = np.array(band_structure);\n",
    "\n",
    "plt.figure(figsize = (5,5))\n",
    "alpha = 0.1;\n",
    "plt.plot(band_structure[:,0], alpha*band_structure[:,1], '.b', markersize = 1);\n",
    "plt.plot(-band_structure[:,0], alpha*band_structure[:,1], '.b', markersize = 1);\n",
    "\n",
    "# plt.plot(Kguesses, alpha*Kguesses**2, '.c', markersize = 0.1);\n",
    "# plt.plot(-Kguesses, alpha*Kguesses**2, '.c', markersize = 0.1);\n",
    "\n",
    "# plt.axvline(np.pi, linestyle = '--')\n",
    "# plt.axvline(-np.pi, linestyle = '--')\n",
    "plt.xlabel('qa', fontsize = 16);\n",
    "plt.ylabel('Energy', fontsize = 16)\n",
    "plt.xlim((-np.pi, np.pi))\n",
    "plt.savefig('Konig_Penny_bands.png', dpi = 300)\n",
    "plt.show();\n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## wave function solutions\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 136,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAD8CAYAAABzTgP2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJztnX2wXlV56H/Pec8HdJwWOMTyGYFbauXKjJYMeupUokTwYwZiRRs6NFGREIS2XKdXk7G0GXoNYP8grVBNkEDS3hGqvWJ61aEEOLUz58UaBwSxgwSqEImSBrDjdTjknPPcP9benH3evB97v/tr7b2f38ye933357PWet79rPWsZ60lqophGIZhhIyULYBhGIbhF2YYDMMwjCWYYTAMwzCWYIbBMAzDWIIZBsMwDGMJZhgMwzCMJZhhMAzDMJZghsEwDMNYghkGwzAMYwmjZQswDMcff7yedtppZYthGIZRKb773e/+p6ouG3ReJQ3Daaedxt69e8sWwzAMo1KIyI/jnGeuJMMwDGMJZhgMwzCMJZhhMAzDMJZghsEwDMNYghkGwzAMYwmZGAYR2SEiz4vI93scFxH5GxHZJyKPishvR46tE5Eng21dFvIYhmEYw5NVuOqdwC3Arh7H3wOcGWxvAT4PvEVEjgP+AlgBKPBdEdmtqi9mJNcS2m2YnoaVK2FqKo8nJKfdhl1Brq1d65dc09MwOQmHDvmVZ77y9596jH/5P4c47/cmueyms8sWx398VX5fKfIFpqqZbMBpwPd7HNsGXBr5/QRwInApsK3Xeb22c845R5MyM6N69NGqrZbq+Ljqhg1uX1nMzDgZxsZUwW1jY+XLFcp29NGqIyNOrpER97tsuXxmwwX7dJRZHeGwHs3/0396+02WYf2YmXF/xFD5JyYsv/qxbZt7QaT8MwJ7Ncb7vKg+hpOBZyO/9wf7eu0/AhFZLyJ7RWTvwYMHEwswPQ2vvALz8+5z2zY4/3xnhIum3XbP3rYNDh9e3H/4sNu3ciVcdVU5soGrxL38MiwsuN8LC+73rl7twQbTbsP73w/b//k05hhjgVFmGefRb71UnoL5TrsNmzcvVf7ZWbfP8utI2m24+mqXXwsLLq+mp3N9ZFGGQbrs0z77j9ypul1VV6jqimXLBo7oPoKVK2F8HETC+5X3sgtfvNolparlGq52G3bsOFI2VbjjDvvfRtm+HX73d+Gee2CBEZw6Ky0WeAcPmjXtRlgr2rPnSCXbs8eMaSehEZ2fX9zXarkXWo4UZRj2A6dGfp8CPNdnf+ZMTcH998OVV8LYmNtXxsuu88U7NgarV7ttYqJ8wzU9vaiDInDWWYsyzc3lXlGpDO02fPzj0f+rMwojzPPX8sdM8ZBZ026ETfeFBRgZgXPPddvIiNv3yiumZCGdRlQERkfhllty72MoyjDsBtYG0UlvBX6uqgeAe4ELRORYETkWuCDYlwtTU/D5z8Pll5fzsguN/9yc+y3iZPnqV9324IPlGq52G555xuleqwVHHQV/8ifus9Vy2zPP2HsOnM6ErraQU04RPr9tlKuuFLOm3ehUsIkJ2LrVbRMTzjiIuIgH40gj+q53wbe+BevX5//sOB0RgzbgS8AB4DCuFXA5sAHYEBwX4FbgKeAxYEXk2o8C+4LtI3GeN0znc5RoR3RRnapJOnQ3bFAVcee1WqpbthQnX7fO+bCjfGKi2DzzlV6BA6/mSRkK5juDoj8y6lytFTnoETE7nzMJV1XVSwccV+DqHsd2ADuykCMuoVspDMcMK3R5ts46jf+qVa710O2Za9fCzp3u/GgtvQj5QtfI8uWLz5uacsfn5hY776enmxldGLbuX3nFVW7f8AZ4/evhk5+M5EcZCuY7/RQMXEz0wsJSd1KT8ytkXTC0q+hw3jjWw7ctbYshpMiKXdIKUdG19EF5YZVgx5YtLg/CloJIn/ywTFvEFCwZOeUHnoWreklnCGteruB2G6691j1nZMS5VAcZ/6kpV6nqrKXnRVjJ/cu/dJ+d8g063hQmJxdd4bAYRda1bIpSsCpgCpaMknWn0YYhDGFttdxnXhFgUTeSqms1J5GviD65OIMqp6Zg0yb3/YYbmtcJ3Wngx8YG6E5RClYFTMGSUbbuxGlW+LZl5UpSdS20LVvKddP0o4g+uSTyNbnFH3UjtVrO1TdQd4pQMN8xBRuOHHSHIjufq8zUVL6t1rCitHXrcHMOFdEn163V2usZSc6tG2El7pVX3Ges/sC8FawKmIINR4m603jDAPnNTRWNYBkfH8512vkyyqNFmeQZRcjjI1ED//DDQ17c1JkITcGS4YG+NN4wZPHy7kUWlZ8iIh+jzxiki0nOrQtRHWm1XH/P3JwLKR6oL3kqWFUwBYuPJ/rSeMOQZ8s1q8pPKE+e+pKk1do070hUR8LRztFopL55Ya4RhylYPDzRl0ZHJUH+nf/r1sEVV6R/kecVvdZuDxcAMux1VSSqI2NjCfWl7OiSsjEFS4Yn+tL4FkNeLdfOFuHatenul4frddhWqyet3ULoFjwACfSlya4RU7DkeKIvjTcMkE/LNesWYagvWc62OqyMnrR2c6ff+ylRepvqGjEFGw4P9MUMQ4QsgwHyCq4I51CK1fGZk4xNCRzp5b4bWkc8iDYpFFOw5HiiI2YYArJuvebRIsyrFZJURk9au7nT+X6anEyhI010j5iCJcMjHTHDEJDlSzdq9MMR/lmQR0Vq2FarB63d3Ol8P6XSkaa6R0zB4uORjphhCMjqpZun0W9qRaosurXqh9aRJrtHjHh4pCOZGAYReTfw10AL+KKq3thx/GbgHcHPXwFeq6rHBMfmcYv3ADyjqhdlIVNSsnrp5m30w4pUGM2XRtYs3JmeuEQzp5eBH1pHmmjVTcGS4ZOOxJlQqd+GMwZPAWcA48D3gLP6nP9HwI7I718kfWaWk+hlTRFzgGXxDF/u4SudE+YVsYperTAF8xIKXI/hXGCfqj6tqq8AdwEX9zn/UtxSoLUkNPp5TiufxWA3X+7hK56MM6oupmCVJgtX0snAs5Hf+4G3dDtRRF4HnA48ENl9lIjsBeaAG1X1ngxkKo0iWr5ZuCJ9uYev+NSqrySmYJVGXOsixQ1EPghcqKofC37/IXCuqv5Rl3M/BZwSPSYiJ6nqcyJyBs5gnK+qT3W5dj2wHmD58uXn/PjHP04ldx4UGW3mi/u2SS5gIyGmYN4hIt9V1RUDz8vAMEwBm1X1wuD3JgBVvaHLuQ8DV6vqTI973Qn8X1X9Sr9nrlixQvfu3ZtK7l6k0cMbboDrrnMt31bLuZOyDFc1iqGbDmT2frIXndGNgvQirmHIwpX0HeBMETkd+AmwBviDLgK9HjgWaEf2HQv8UlVnReR44G3AZzOQaSjS1viLbvnaOyZ7uukAZNQS9GgAk+ERHupFasOgqnMicg1wLy5CaYeqPi4i1+N6wHcHp14K3KVLmyhvALaJyAJuptcbVfUHaWUalrShpkX6pT3UpVrQq78zkxBkjwYwGR7hoV5kMo5BVb8BfKNj3593/N7c5boZ4OwsZMiCLGr8RQ3Y9FCXakEvHcikJWidqUY3PNQLG/kcIU2Nv2i3ThpdylrWOrm0eulAJi3BpoQ6mYIlw0e9iDPYwbfNtwFuZY3DmZlxA6+SPC9rWW0MkrEEUzCvocABbo2nrHE4U1Mu6ilJBSNrWes2BqmpC4dlhilYLTBXUgaU7SJM0tLOWtay054l1qGfAaZgtcAMQwaU6SJM+jLLWlYf3aPDYh36GWAKVgtSD3ArgzwHuIVUpb/LBtVlR68xDLnoQVUUzMifAnWhyAFutSNpLbzM/7i1tLOjs3IKObmWzGdlhHiqC2YYupDEpVB2uVpLO1ui41BuuCEn15L5rIwQT3XBDEMXktTCfSjXJq6CmDXdWn25tcasmWeEeKoLZhi6kKQW7mm5Ggno1erLrTVmzTwjxFNdsM7nDPChHzGuDHnK6kM+DIN14GeIKZjXWOdzgZTtyonbz5Fnf0jZfS1psFZfRpiC1QYb+VwD4g4OzXMQaZUHqBaxHGsjMAWrDdZiSIEvLdu4Nd48a8ZVr3WX3eqrBaZgtcH6GIbEt5at9TEYXmAK5jWFLe1ZBkUZhn56aB2W9aFXORfyHrKXXXMpoezjGoZMpsEG3g08AewDNnY5/mHgIPBIsH0scmwd8GSwrYvzvCKm3R4026/NBlwPepVjIeVrStRcSip7ipp2W0RawK3Ae4CzgEtF5Kwup96tqm8Kti8G1x4H/AXwFuBc4C+CdaBLZ1Bfl68dljZtdDJ6lXMhfZ3WodpcPC/7LDqfzwX2qerTACJyF3AxEGft5guB+1T1heDa+3Ctjy9lIFcq4vR1+dZh6Vu/RxXoVc6F9HVah2pz8bzsszAMJwPPRn7vx7UAOvmAiLwd+CHwP1T12R7XnpyBTKnxdEBiX3yYnqNq9CrnQsq/ikpmZIPnZZ+FYZAu+zp7tP8J+JKqzorIBmAn8M6Y17qHiKwH1gMsX758eGkT0KtF4Gt/Yb9KSFEy+5o3/ehVzoW0CH1rdg7LMAXfbsOuXe772rXxrquigvXC47LPwjDsB06N/D4FeC56gqoeivy8Dbgpcu3Kjmunuz1EVbcD28FFJaUROA0+u2t6VUKKktnnvDFyJMnQ++lpmJyEhx+G22+Hw4fdsdtvh/e9D044obeRMAUrjCwMw3eAM0XkdOAnwBrgD6IniMiJqnog+HkR8O/B93uBLZEO5wsAr4M+fXfXdKuEFCWz73nTSZ0qn6USp+C3b4drroG5OVAFEfcZcvgw3HOP+3777XD55UcaiKopWIVJbRhUdU5ErsG95FvADlV9XESux4VG7Qb+WEQuAuaAF3Dhq6jqCyLylzjjAnB92BHtK573GXWlKJmrlDdW+cyQQQXfbsPVVzujENJv/NThw7BtG+zcubRgqqRgFccGuA2B7zXNbvJZH8NSbIBixvQbJbh5M9x331JjMDICo6Pw3ve631//+qJbKUQErrwSPv/5wc8xYmEjnxuK1YTjYflUAN3cR60WfOITcMwxR9Zcdu2Cn/50qZEYG+vuVjKGwqbdzpAqVVLMDRuPftGChZd3lRQsLp3uIxF417tc62FQCNhVVzlXkmpvt1LV8b3M4wyP9m0rYkqMkKrNWlA1eX2j8PyrY4HNzKhecIGqiKp7vauOjcVPW5gn0etFVDdsyFfuoiixzClqSoy601kD37XL7yknfJ2qoyoUPlOB51MjJCb00e3Zs+g+Gh2FW26Jr4yhEl95pXMlgbvXHXf4+8dLQgXK3FxJA4gGQoyOwo4drjx99kt7PG7GewoPfKlbpM2uXfDyy+5FPjICq1b1dh/1I6rEUbfS5s3D3c8nKlDm1vkcg9Ad+MwzcNttFslSB/q5eK2PYUjabZeGV15xvycm4MEH06UpbIHMzsLCgjM2ExP+1sriUlKZFzrtdtFbkX0MUeroDm4iVo45sWWLy9Ss+wTCPouREXfvVss9y0gMMfsYzJWUAM/nvTqC6AwEhw4VK7PPlWCL3MqYaKjpaPBKGR93IaZZMDXl3Ef/+q+u5SACL73kOvt8VLA6EMd6+LaV1WKoEmGtOKxkjYwUVzv2vUbuu3yVYmZGdWJiMXpodNS1FPLI1G3bXHRTGK1UpFLXBCwqqdmEteKFBfd7YaG4AAjfgy4scitDwsIOmZ+H5cvzydRDh5wih/2iRSp1wzBXUk0JAx+ifXZFBUBUIOjCIreyYuVK5z6KjlTOezKuUKnDUFgfFazimGGIic8+825E+0OK7mOoWl+MkRIRt42MwOc+l1+Bh4q1a5ebgTWcasPIHDMMMajqvDpl1oqtRt4Qpqed+yh8QR861Pf01ExNuWeGLqX5eYseyAHrY4jB9LRrvc7Pu09zaVaXdrv/yPVBx3Ol1IcPQbvtBveMjrqBPUX7KkdGXEtlcjL/Z2ZFRcrYWgwxmJxc2olbJT00FhnU8iu1ZVi1ZmlU3lYLrriiuBlQp6Zg61Y3c+v8PFx7LZx9tt/5BZUqY2sxxODQIVc5AfeZd2s5aypSScmdQdFSpUZT+R7K1UlU3jwjkXoRRihVKTKpQmWciWEQkXeLyBMisk9ENnY5/gkR+YGIPCoi94vI6yLH5kXkkWDbnYU8WbNypRuF32q5zyoFQYSVlOuuc59NNg6hB6KX12PQ8VKF843JSVdLKjLcLUo0v1ot59LyXbmrVMZxBjv023DLeT4FnAGMA98Dzuo45x3ArwTfrwLujhz7RdJnljHAbWbGjcKv2lia6CwFNpPA4HIstZyromTR0ZOjo27gWVlybNjgBthVZbRiyWVMgVNinAvsU9WnAUTkLuBi4AcR4/Ng5PyHgMsyeG6hVDXKpgpjCopkUDmWWs5VUbLo6EmR8nyrYYTS3Fx15jepSBln4Uo6GXg28nt/sK8XlwPfjPw+SkT2ishDIrK610Uisj44b+/BgwfTSZyQqvrow7EXW7eWM8q3qvlm9CEaiVR2VFAoS+hOqopLqQrEaVb024APAl+M/P5D4HM9zr0M12KYiOw7Kfg8A/gR8N8GPdNWcBtM2XKX/XwjB6KFOjrqPsuarygqy8SE6urV1XIplQQFzpW0Hzg18vsU4LnOk0RkFfBp4CJVnY0YpueCz6eBaeDNGciUGRUKJFhC2XKX/fxuWAsmJZ2RSGVGBUVlmZuDX/7ySJeSMTRZGIbvAGeKyOkiMg6sAZZEF4nIm4FtOKPwfGT/sSIyEXw/Hngbkb4JH6hSIEGUsuUu+/mdWHRWBkQLdWzMLwX7wAf8UriKk7rzWVXnROQa4F5chNIOVX1cRK7HNVt2A38FvAb4sogAPKOqFwFvALaJyALOSN2oql4ZhqrO+9MpNxQ7fb1v+WZrMGTEunXuM1xroawC7qZgZ5/t5lEyUmNLezaACg24zA3Lg5RUIQOrIGPJxF3a00Y+J6SKfmof/f1FM2gNBm/K1RtBOqiCEvkso6/l2gObKykBVa2Q2FgGR68Qcm/K1RtBusgVhqiCv0oUVfRo6GrZeehrufbBWgx96DTyPldI+mErlvXHm3L1RpAI4UvtttvcNNdXXOGvEoWKfsUVbnzFbbf5EWngY7kOwFoMPehm5Ktc867IgMtS8KZcvREkQvSlBsVPlpcUH0dD+1iuAzDD0INuRn7TJr8ibYx4DFp9z5sIKm8EiVDBl5p3MvtYrgOwqKQeVNAtaHTByjEFoUUtem3YLGi3F0NXi1onogLEjUqyFkMPKmjkX6VfDbmMtavLXC/bxi8MSZxVjXz5c/SSZedOJ//OnVYjSIgZhj5U0S/f7/9cRu257Bq7b16FytDPopZdqFF6yWI1glRYVFLN6BcAUUZwRNkBGRaRNST95jQpu1Cj9JLFtzlZKoa1GGpGvxpyGbVnH2rsVWz5eUF0+otoBvpQqINkifqCJycXDYYpQiys87kHPrlQk2J9DEYq4riKfCrUQQrvi9vLA6zzOQVV16V+NeQyas9WY68YcfzzPhVqP1msr2EorI+hC4NcqBWb9qTRxCkrL8uzTKHq5J/3IS1eKlh/rMXQhX4u1Kq3JppEXI+Id+VZplDR9WCrNnahG2FfQ1nTcXupYIOxFkMX+kWy+BSQMSwVrMAMRZyy8rI8yxIquprRtddW3yhE2bmznLmTvFSwwWRiGETk3SLyhIjsE5GNXY5PiMjdwfFvi8hpkWObgv1PiMiFWciTBVNTbgqMzv+FDy3TNDRpJbM4ZeVleZYlVEVfYgMpM11eKthgUruSRKQF3Aq8C7f+83dEZHfHSmyXAy+q6m+IyBrgJuD3ReQs3FKg/x04CdgjIr+pqvNp5cqLKo+Ihmb1xcUpKy/LsyyhfApDzZIy0+Wlgg0mdbiqiEwBm1X1wuD3JgBVvSFyzr3BOW0RGQV+CiwDNkbPjZ7X75m2gtvwVNTlaRRBnecXqvK8TxlSZLjqycCzkd/7gbf0OidYI/rnwGSw/6GOa0/OQCajBxWtwBh501ljCNd0rguholutKBZZ9DFIl32dzZBe58S51t1AZL2I7BWRvQcPHkwoYnyq3jEbR/5e/Sdly2WUSFw/vM8FOUi2uvah5EAWLYb9wKmR36cAz/U4Z3/gSvo14IWY1wKgqtuB7eBcSRnIfQRVd7P4Kr+vchkR4vjhfS7IOLLVtQ8lB7JoMXwHOFNETheRcVxn8u6Oc3YDwcQrXAI8oK5zYzewJohaOh04E/i3DGQaiqpXKHyVvwy5fK7Yekd07EK/2QZ9VTCIJ1t06c9wHiijK6lbDEGfwTXAvUAL2KGqj4vI9cBeVd0N3A78nYjsw7UU1gTXPi4i/wD8AJgDri4zIqnqFQpf5S9aLp8rtt6RJLN8VTBIJput0zCQTEY+q+o3gG907PvzyPeXgQ/2uPYzwGeykCMtVe+YTSp/UfOgFZ2vTQrJTU2SzPL5DxJXNlOOWNjsqg2lzrXqOqctc5qWWU1Lbwc2u2oB+DTzcFLqXHGKW3n0vvyKErDXugt1pMh1GrxXsD6oauW2c845R8tmZkb16KNVWy33OTNTtkTJqLr8afE+/UUI6H0m5Ejeafc0b3H9vgPfsTaJXoQkkSw+B2jEoelLXnpffkUI6H0m5Ejeaa943porKSCp69HnAI24+LTWStF4X35FCOh9JuRI3mmveN6aYQhI6nP3OUCjycR163pffnkLWLd1F5KS9zoN3itYfywqKaDhwQq1wMowJpZRjgbmQ9yoJOtjCKi6z33Ykb5FjxDO83kVd+sWx7AZVYXh5E3qKMwRcyVFqKrPfdiKT9EVpryfV3G3bnFMTsLICKjGz6gq1K7TdBS2WvDMM+4evqWrBKzFUAOGrfgUXWHK+3lVb/UVQrvtlu2cn3fGYevWeBlVhdp1UhmjcyeJlLP0p6dYi6EGDFtTLrqGXcTzqtrqK4zw5bmw4F6Ghw7Fu64KzbFhZJyacnkyN1fP0Z5DYoaBag9QhOEDIIoOnKh4oEY9GPYFX4XCG1bGKhi9gml8VFIVXKeGkQm2vGVv6rysaQSbKykmWcwZVPUWR11IUg6VKrMshLUa0GCynI67Ugp2JI03DGlbkfZ/84Mk5VCpMstK2DrPmpgFWeZPpRSsO42PSkobyVKFYI04VCFEvR9JyqFSZZaVsGENqNUyP3o3ssyfSilYd1K1GETkOOBu4DTgR8CHVPXFjnPeBHwe+FVgHviMqt4dHLsTOA/4eXD6h1X1kTQyDUOaSJY69FvVoIKTqBwqVWZZCtuk6bWTkuV03HkpWJHuqThTsPbagM8CG4PvG4Gbupzzm8CZwfeTgAPAMcHvO4FLkj7Xh2m3o8zMqG7Z4s3MuonZssXNDgzuc8uWsiUajiTlkEuZ5aUIae/r6RTQXpJVXmWtCxnJRcxpt9P2MVwMrAy+7wSmgU91GJ4fRr4/JyLPA8uAl1I+2xuqHjtfqRp0H5KUQ+ZllmezK62w1r8Qn6zyKmsFK7gM0/Yx/LqqHgAIPl/b72QRORcYB56K7P6MiDwqIjeLyERKeRJTdd86pE9D0SOG65DnR+CrX7nddlM9jI4O7z+vYoENK3O0ryE6TUbZFN1HNKhJAewBvt9luxh4qePcF/vc50TgCeCtHfsEmMC1OP68z/Xrgb3A3uXLlw/VjOqkDi3sqqWhavLGJpqw8XHVDRvKT1wWMlWxwNLKPDPj8mpiwp90hzKl1CuyWsFNVVep6hu7bF8DfiYiJwIEn893u4eI/CrwdeDPVPWhyL0PBPLOAncA5/aRY7uqrlDVFcuWLRskdix8reQloWppyENeLyq0Ps67E83s+XlYvrwZYXdpZZ6acnnVOU1GWYRuyttuc2MsCiCtK2k3EIQ6sA74WucJIjIOfBXYpapf7jgWGhUBVuNaIoVRhwi+qqUha3nD/8x115X/HvbuhRLOojoyMnxmV03BIBuZw3uMjDhDPzmZsZAJKMM4x2lW9NqASeB+4Mng87hg/wrgi8H3y4DDwCOR7U3BsQeAx3AG4e+B18R5bpZRSVWPKFLNLg1F5UWWz/Euoip0Y4yMqI6Oqm7bVn05qvgnyULmbdtUx8ZcHpblTsrYrUVMV1Iqw1DW5lu4ah2ooitZ1VO5fXiheGcxK0jZeZhDv1Vcw9D4kc9Z4YWfOwVVdCVD8oiqQsrp0CE3rfXCQjkdKVlEIhnDRShlqWBZ9BENSxzr4dvmW4vBy1prQuqQhkEUlsasa3pJBPcxOqrKJHHlZK1gObQ8sRZDf/Iy7FWqbUdpwupnhZVT1hFKw04EVXQts44kCSjIUsGGXWkvIxo5u2rWg1SbOHK4ihRaTlmuDJZE8GHWczb6E+b/7Gz/CKUsFSw0MklX2suIRhqGrEeXV2Fxq7qSZF6xwssp7gtlEHEFL7mWWVumplxeXnONy9trr4Wzzz4yb7NSsGgfEZRj4OP4m3zb0vYxNMGf3gQqUY5FRSjNzKhecIF7jkUiZU9REUo59xFhfQy9qYs/Pa8Im6IirNI+pxJ9O9EIpdlZ2Lw5nwI7/3zYs8c9J82Ats77VjnULkv5i5pDadcuePnl8vuI4lgP3zbfopLKIK/aclG18CyeU4kWQ3SgGeTTcojWZkdGXMshi9GO3mduH/KQP+85lGZmXCvB9RC552Sc71iLod7kVVsuqhaexXMq0fILhVy1ytXkFxZcjTBceD4tnWMWJiZcqyRtZlSiOdaHPOTvjFDKuhw3b3b3Btcn9ZGPlKfUcayHb5u1GKzFUDnyqA3m6Y+uegHl+QfJqxzzbFUGYFNi9KaKU790w9cFw3x7jjds2KAqkp3LJ3q/PDpEq15AecmfZzlm5QrsgRmGHlS9ImRUmKxqhqGve2ws25qrEY/OchQZfqLCAvoVosQ1DI3rY8jbdVr1QI6qMEw+l1423fobkkYqbd8Ob387fOELcPiw21e2P7ppRMtRxL3S5+bg4x+Hq66KX5a+9StEiWM9fNt8bTFYa6QYhslnr8pm2BrnzMxi9FG4iXiQoIYyM+PKLVoecWv94fiWqAupgHLEWgzdyTOSpeqBHFFKr133YZh89qpsktY42223/2MfcwkIGRmBK6/0OCSr5kxNwa23wtjY0v2zs7BxY/e7XSGhAAAS6ElEQVQ/UFiWH/+4a/GpOh1YtcqvcoxjPXzbfI1K8qpWmgLf01H5FkNIrxrn6KjqJz/p+hFWr17alxDdVq8uOwWGqivHc889snzClmC0LCcmFlsJ4TY2VphCUkTnM3AccB9uBbf7gGN7nDfP4uptuyP7Twe+HVx/NzAe57m+GgbV6gdyqJa/PkkchslnL8um06UQdxsf9ywhDaebmy/OVvAqf0UZhs8CG4PvG4Gbepz3ix77/wFYE3z/AnBVnOf6bBjqgJe16zoTRhnFebGMjdkaC74S18iLlLZWRlzDIO7c4RCRJ4CVqnpARE4EplX19V3O+4WqvqZjnwAHgRNUdU5EpoDNqnrhoOeuWLFC9+7dm1jeJDNx+k7eaSkqr+pUJqnZvt3N4Dk3514hIWNj8L73wQknwNq1xWRUnQqmyLSEz3rpJbj55u5lefnlxZVjByLyXVVdMfC8lIbhJVU9JvL7RVU9tst5czg30hxwo6reIyLHAw+p6m8E55wKfFNV3zjoucMYhqzXYCiTuqRl2HTU6Z11BGHiJifh4YfdvqJfInVRMCg3LT6UZQdxDcPA9RhEZA9wQpdDn04gz3JVfU5EzgAeEJHHgP/qcl5PKyUi64H1AMuXL0/waEfWazCUSV3SMkw66vTO6ooPqyXVRcGg3LT4UJZDMjBcVVVXqeobu2xfA34WuJAIPp/vcY/ngs+ngWngzcB/AseISGicTgGe6yPHdlVdoaorli1bliCJjuisuVVf2KrItOQZtjpMOrwKO60r9mdpPGlXcNsNrANuDD6/1nmCiBwL/FJVZwP30duAz6qqisiDwCXAXb2uz4o6rbJWVFryrp0Pk466LKPqNfZnaTxp+xgmcZFFy4FngA+q6gsisgLYoKofE5HfAbYBC7gWylZVvT24/gycUTgOeBi4TFVnBz132M7noqiLD/yGG+C661ztvNVygwI3bSpbquT5W5fyeJXaJajiVKg84vYxpApXLWvzOVy1TqGedUhLHdKwhNolqOJUrDywKTHKoU4+8EoshDOAOpUHUMMEVZyalkfaPgajg7r5wH0KrBimxV638qhfgipOTcsjVR9DWVgfQ/NI0xFeu/KoXYIqToXKI7NxDEZy8qxll6GDRTxz0DPShKP71OrJhLopWBHkma7aKZgZhkpRxuCuIp4Z5xk1bbH7RV1HD9Y1XTlinc8Voox+riKeGecZdegI956adqTWNl05Yi2GClFGrbmIZ8Z9Rg1b7H5R12ZZXdOVI9b5XDGa2sdgFERdC6Ku6UpIIbOrlkVVDIPpYnYMm5e1LoNaJ64iVKwMLCqpZOrY31XWfyDN9Nx1K4NXqXXiKkKNy8A6n3Oibv1d4X/guuvcZx6zrfZi2LysWxksodaJqwg1LgMzDDlRt9l+y/wPDJuXdSuDJdQ6cRWhxmVgfQw5krXrpUx3ZlGt5l5ptD6GLtRJwYqi4Xlmnc81wwd3ZhHrTJedxsbShMxvQhoHENcwmCupIvjgzpyacusx5PVf6pbGPFeQMyL4oGB504Q0ZoRFJVUEn8bo5NVy6Ezj5GTjK3jF4ZOC5UUT0pgRqQyDiBwH3A2cBvwI+JCqvthxzjuAmyO7fgtYo6r3iMidwHnAz4NjH1bVR9LIVFd8WaEwz9Z4ZxrTTJxnJMQXBcuTJqQxI9K2GDYC96vqjSKyMfj9qegJqvog8CZ41ZDsA/45csr/VNWvpJTDW7KsXfswJUTeL+vONA5bwatYn+Dw1E3B8iarNNZcwdIahouBlcH3ncA0HYahg0uAb6rqL1M+txLUsa+ryNb4sBW8OuZ7VxqTUM9oQL6n7Xz+dVU9ABB8vnbA+WuAL3Xs+4yIPCoiN4vIRK8LRWS9iOwVkb0HDx5MJ3VB1LGvq8hZToetlNUx37vSmIR6RgPyfWCLQUT2ACd0OfTpJA8SkROBs4F7I7s3AT8FxoHtuNbG9d2uV9XtwTmsWLGiEjG2de3rKsLjkKZSVtd8P4LGJNQzGpDvAw2Dqq7qdUxEfiYiJ6rqgeDF/3yfW30I+KqqHo7c+0DwdVZE7gD+NKbclSCrvi4f3Zl5ytRuw+bNMDsLCwvDrdjWiD7GOitYnqRNbxMUTFWH3oC/AjYG3zcCn+1z7kPAOzr2nRh8CrAVuDHOc8855xxtCjMzqkcfrdpquc+ZmbIlylem8N4jI6rgPn1Jdy3xUcHypGnp7QDYqzHesWn7GG4E3iUiTwLvCn4jIitE5IvhSSJyGnAq8C8d1/9vEXkMeAw4HvhfKeWpHT66M/OUKbz3wgKMjMCqVbXs2/MHHxUsT5qW3iFJFZWkqoeA87vs3wt8LPL7R8DJXc57Z5rnV4lhW68+ujPzlGnlShgddYZhbMy5lGxupBjUScHyJIv0NkDJbORzAaTpSPXRnZm3TOH0XUmn8WpAFGF36qZgeZI2vQ1RMjMMBZB2UJiP445CmcK5jLJ6p0xPu3xSdZ9J8qqxI6XrqGB5kia9DVEyMwwFUNfWeh6VpzR5Vdd8HkhjE14CDclrMwwFkKb16rM7M4/KUzSvJicX+wbj3LdpXpFXqauC5c0waW+Iktl6DB7juzszD/nC/+rkJFx7rb9prwW+K1ieNDTtcddjsBaDx/juzgwrT7t2ZXO/6H9VxEUmDTPAzYiJ7wqWJ01OewxsoZ4CSbroTFWWlN25E267zb3U0yyoE/2vLiy4dCdJuy3qQ7JMqIqC5cEwaW+QglmLoSCGablWwZ2ZZcWrs19v61Y4dChe2hvqGVhK0kyogoLlRdK0N0zBzDAUxLAvUN8jCbMK0gj7FpIYgyjmGWC4TPBdwfIkSdobpmBmGAoi6Qu0KsEiaaKIQvpVxuLmQ0OiCPuTJBOqomB5YwrWnTgTKvm2VXUSvZkZ1S1bBs/bVcV5vtLIvGWLuw7c55Ytw90zbv7WmjiZUEUFy4MGKhgFTaJnJGBqCjZtct/79WFVcZ6vYWVut+GZZ9z8SJ39gHHvGfYJgsvfJleAmZpaXDC7TgqWB0nyoWEtLHMlFUycPqwqtlpDmWdnXajp5OTga6J50WrBFVfA2rWL+REnHxrWJziYuipYHsTNhwYqmRmGgonbh7VunfuMvih9ZmrKdRxfc41L27XXwtln95d91y54+eXFyfKWL196fpz+i4b1CQ6mrgqWB3E7yBqoZGYYCmZQzbqzcrJ2bSliDsWhQ4uD0mZn3ZTZvabNbrdhx45FozA62r3CFl7bq8I2OenWbVBtduX3VaK14FbL+ena7aW9+VVVsDwYpGBRXyc0RslS9TGIyAdF5HERWRCRnsOsReTdIvKEiOwTkY2R/aeLyLdF5EkRuVtExtPIUwXCmnWr5V6g11671BVcZfdv+E4aGXFp27On+6C3cOnOuTn3WwQ+8pHelbBeedJuu/ybn3fP3Lq19hW5wYS14CuucBnbOfKwygqWF/0U7PzzXR6qujxtgBsJ0o98/j7we8C3ep0gIi3gVuA9wFnApSJyVnD4JuBmVT0TeBG4PKU8laBbzbrd7t8RWwXCd9KqVYvG4eWXl06ZsX07nHce3Hff4iptRx3Vv+IaNThhK6tzXWhVl68GriCWL3eWd35+sRCqrmB50U3BYNHXOT/vtk5fZ52JE7o0aAOmgRU9jk0B90Z+bwo2Af4TGO12Xr+tquGqIZ3rGou4iLmxMfc5Pq66YUN1o+JmZlwa3OvapWv1areFYalhui+4IF46t21z9xkZUR0ddfcRsXWhe9JZCK2W20ZGqq9geRBVsPFxp6xjY4v5NzFRi/zCo3DVk4FnI7/3B/smgZdUda5jf+2J1qxFFhelOXy4HpWTqSn46Edd2sCl65573DY/v3je6Gj8pTujraywIqzqnmHrQnehsxBCxQozsMoKlgdRBXvlFaeshw+7Y4N8nTVkoGEQkT0i8v0u28UxnyFd9mmf/b3kWC8ie0Vk78GDB2M+2l+mptxLsdVaul+kHi38tWudi0i6lTLOKNxyS/z/Wtja77xfEuPSOHoVQqtVfQXLml4KJjLY11lDBhoGVV2lqm/ssn0t5jP2A6dGfp8CPIdzIx0jIqMd+3vJsV1VV6jqimXLlsV8tN9MTcGtt7pF70dGnF5eeWU9ar9hq+jKK136QsbGYMMG+Na3YP364e43MeHyK6lxaRyWafHpzKuwD6Yuf8iEZLJQj4hMA3+qqkesnhO8+H8InA/8BPgO8Aeq+riIfBn4R1W9S0S+ADyqqn876Hl1W6in7oMq2+3FDugswubrnl+5YJkWnxrnVdyFelIZBhF5P/A5YBnwEvCIql4oIicBX1TV9wbnvRfYCrSAHar6mWD/GcBdwHHAw8Blqjo76Ll1MwyGYRhFUIhhKAszDIZhGMmJaxhsEj3DMAxjCWYYDMMwjCWYYTAMwzCWYIbBMAzDWIIZBsMwDGMJlYxKEpGDwI+HvPx43OA63zC5kmFyJcPkSkZd5Xqdqg4cIVxJw5AGEdkbJ1yraEyuZJhcyTC5ktF0ucyVZBiGYSzBDINhGIaxhCYahu1lC9ADkysZJlcyTK5kNFquxvUxGIZhGP1pYovBMAzD6EMtDYOIfFBEHheRBRHp2YMvIu8WkSdEZJ+IbIzsP11Evi0iT4rI3SIynpFcx4nIfcF97xORY7uc8w4ReSSyvSwiq4Njd4rIf0SOvakouYLz5iPP3h3ZX2Z+vUlE2kF5Pyoivx85lml+9dKXyPGJIP37gvw4LXJsU7D/CRG5MI0cQ8j1CRH5QZA/94vI6yLHupZpQXJ9WEQORp7/scixdUG5Pyki6wqW6+aITD8UkZcix3LJLxHZISLPi8j3exwXEfmbQOZHReS3I8eyz6s4639WbQPeALye/mtRt4CngDOAceB7wFnBsX8A1gTfvwBclZFcnwU2Bt83AjcNOP844AXgV4LfdwKX5JBfseQCftFjf2n5BfwmcGbw/STgAHBM1vnVT18i53wc+ELwfQ1wd/D9rOD8CeD04D6tAuV6R0SHrgrl6lemBcn1YeCWLtceBzwdfB4bfD+2KLk6zv8j3FIBeefX24HfBr7f4/h7gW/iVr58K/DtPPOqli0GVf13VX1iwGnnAvtU9WlVfQW3LsTFIiLAO4GvBOftBFZnJNrFwf3i3vcS4Juq+suMnt+LpHK9Stn5pao/VNUng+/PAc/j1gfJmq760kferwDnB/lzMXCXqs6q6n8A+4L7FSKXqj4Y0aGHcKsl5k2c/OrFhcB9qvqCqr4I3Ae8uyS5LgW+lNGze6Kq38JVAntxMbBLHQ/hVr88kZzyqpaGISYnA89Gfu8P9k0CL6nqXMf+LPh1VT0AEHy+dsD5azhSKT8TNCVvFpGJguU6Sty62w+F7i08yi8RORdXC3wqsjur/OqlL13PCfLj57j8iXNtnnJFuRxX8wzpVqZFyvWBoHy+IiLhEsBe5FfgcjsdeCCyO6/8GkQvuXPJq9HBp/iJiOwBTuhy6NMabz3qbsvUa5/9qeWKe4/gPicCZwP3RnZvAn6Ke/ltBz4FXF+gXMtV9TlxK+89ICKPAf/V5byy8uvvgHWquhDsHjq/uj2iy77OdOaiUwOIfW8RuQxYAZwX2X1EmarqU92uz0GufwK+pKqzIrIB19p6Z8xr85QrZA3wFVWdj+zLK78GUahuVdYwqOqqlLfYD5wa+X0K8BxuHpJjRGQ0qPWF+1PLJSI/E5ETVfVA8CJ7vs+tPgR8VVUPR+59IPg6KyJ3AH9apFyBqwZVfVrcOt9vBv6RkvNLRH4V+DrwZ0EzO7z30PnVhV760u2c/eLWOv81nHsgzrV5yoWIrMIZ2/M0snxujzLN4kU3UC5VPRT5eRtwU+TalR3XTmcgUyy5IqwBro7uyDG/BtFL7lzyqsmupO8AZ4qLqBnHKcFudT06D+L8+wDrgDgtkDjsDu4X575H+DaDl2Po118NdI1gyEMuETk2dMWIyPHA24AflJ1fQdl9Fed//XLHsSzzq6u+9JH3EuCBIH92A2vERS2dDpwJ/FsKWRLJJSJvBrYBF6nq85H9Xcu0QLlOjPy8CPj34Pu9wAWBfMcCF7C05ZyrXIFsr8d15rYj+/LMr0HsBtYG0UlvBX4eVHzyyas8etjL3oD34yzpLPAz4N5g/0nANyLnvRf4Ic7ifzqy/wzcH3cf8GVgIiO5JoH7gSeDz+OC/SuAL0bOOw34CTDScf0DwGO4F9zfA68pSi7gd4Jnfy/4vNyH/AIuAw4Dj0S2N+WRX930Beeauij4flSQ/n1BfpwRufbTwXVPAO/JWN8HybUn+B+E+bN7UJkWJNcNwOPB8x8Efity7UeDfNwHfKRIuYLfm4EbO67LLb9wlcADgS7vx/UFbQA2BMcFuDWQ+TEi0ZZ55JWNfDYMwzCW0GRXkmEYhtEFMwyGYRjGEswwGIZhGEsww2AYhmEswQyDYRiGsQQzDIZhGMYSzDAYhmEYSzDDYBiGYSzh/wMcQONctr5SGQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1bd2ddd4a90>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "## do a scan of K...\n",
    "Kguesses = np.linspace(1e-3,4*np.pi, 4);\n",
    "band_structure = [];\n",
    "for kguess in Kguesses:\n",
    "    val = RHS(kguess);\n",
    "    if(abs(val) <1):\n",
    "        q = np.arccos(val);\n",
    "        E = kguess**2;\n",
    "        band_structure.append([q,E]);\n",
    "        x1 = np.linspace(0,1, 100); #a has been 1 in everything we've done\n",
    "        x2 = np.linspace(-1, 0, 100);\n",
    "        C = 1/2; D = 1/2;\n",
    "        A = C*np.exp(1j*q)*np.exp(1j*kguess);\n",
    "        B = D*np.exp(1j*q)*np.exp(-1j*kguess);\n",
    "        psi1 = A*np.exp(1j*kguess*x1) + B*np.exp(-1j*kguess*x1)\n",
    "        psi2 = C*np.exp(1j*kguess*x2) + D*np.exp(-1j*kguess*x2)\n",
    "        plt.plot(x1, psi1, '.r');\n",
    "        plt.plot(x2, psi2, '.b');\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## negative sign of the potential"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAU8AAAFFCAYAAABsVm+UAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xl81fWd7/HXJ4EQEsKigAuIqBXcRlTi3lGnase2VvGOvVO30ar1+ujeaXtttVWx1baP3plxtFPvtCpqa9W61O1OrdQZFXEFi8uIIFoBFSWIQGLI/rl/fPLznOw5hyRnez8fjzxyzu/8cvKFnPM+3/1n7o6IiGSmLNcFEBEpRApPEZEsKDxFRLKg8BQRyYLCU0QkCwpPEZEsKDxFRLKg8BQRyYLCU0QkC6NyXYBsTZ482WfOnJnrYohIkVm6dOkGd58y0HkFG54zZ85kyZIluS6GiBQZM1s9mPPUbBcRyYLCU0QkCwpPEZEsKDxFRLKg8BQRyYLCU0QkCwpPEZEsKDxFRLKg8BQRyYLCU4pSRwfccQe0t8Of/xxf7e1xrKMj16WTYlCwyzNFIILwtttg8WLYaac4VlYGZnDZZfDmm3DttXH8q1+F738fnngCDj00HnOHdevgyCPhtNPiZ0UGQ+EpBae1FU4/PWqS77wDzzzT85yJE+FHP4Jvfxs++ck4tv/+cf4118DPf971/Ouui5DdeWcoL4ff/hZGjx7+f4sULivU67bX1ta6NgYpHW1tUXNcvjy+1q/v+vixx8LRR8ftsjI44QQ46KCogaZLb84nNc/HHoNHHul63tSpMH06nHwyXHJJBKqUBjNb6u61A56n8JR81tISQfjyy1BX1/WxY46ByZNh3rxta3InTf9774UNG+DRR7s+Pm1aNPNvv1210VKg8JSC1twMRxwBr74KjY2p4zNmwKxZcNVVUFvbs2a5rdxhyRK4+GJYuhQ++CD1WE0NzJ0Lf/wjVFQM7e+V/DHY8FT3uOSV5uZobo8fD88/nwrO3XaLPsw33oCFC+Hgg4c+OCGe8+CD43esXw9f/nJqIKq+PmqlEyZErbelZeh/vxQOhafkhZaWCKTx42NaURJMu+8Ot9wCq1aNfN/jqFExsLR2LfzwhzGYBNDUFP2kEybAKafEAJaUnhENTzO70czWm9nLacd+ZmavmtmLZvZ7M5s4kmWS3OrogAULYNy4CKQkNKdOhaeeitA866zcTiEqL48pTmvXws03x0g+RIjee2805598Mpr8UjpG+iV5E3BCt2MLgf3cfX9gJfC9ES6T5MjWrRGS556bqr1VV0dArVsHhx02PE3zbJWVwT/8QwwqXXEFjBkTx5ubY57ozJkRqFIaRjQ83f1xYGO3Yw+7e1vn3aeB6SNZJhl5ra0xBai6Gt5/P3X8wgth06YIqHyerF5eDj/4QfSBnnpq6viaNfFvOv/8mFolxW3ER9vNbCbwoLvv18tjDwB3uPtv+vjZC4ALAGbMmDF39epBXadJ8siHH0azNz1c9tgjpiJVVuauXNti69aodabPPa2sjA+GqqqcFUuyVHCj7WZ2CdAG3NrXOe7+S3evdffaKVMGvDKo5JGWFjjqqOjbTIKzvDz6Cl97rXCDE2Ds2OhmuOmmVDdDU1PUQk8+WQNKxSovwtPMzgZOBM7wQp14Kn1qaIiAWbQodey442Ia0uGH51e/ZrbKyuDss6MWul9am+r++6P2+eGHuSubDI+ch6eZnQBcBJzk7o0DnS+Fo7UVTjopRqOTnYwqKiJMFy4szonmY8bAiy9GjTqZVtXWFjXu885TX2gxGempSrcBTwGzzewtMzsP+DlQAyw0s2Vm9n9HskwyPBobIzAeeCB17NRTIzirq3NXrpFgFjXqxsZY1pm48UaYNClqp1L4Rnq0/TR338ndR7v7dHe/wd0/5u67uPsBnV8XjmSZZGi1t8dIdHV1as5mWVmMTN95Z2mtDa+oiLmqixenjjU0RDP+hhu0r2ihy3mzXYpHUxPssEMso0zsvnuqFlqKzGKN/tatkD7Gef75sM8+WuJZyBSeMiSSQaH0eZsLFsRIejKZvJRVVsaI/Pz5qWMrVmgwqZApPGWbtLXFQEhNTepYVVXUNs85J78nu4+08nK49NL4oBnVuQ15e3vUyq+/Xs34QqOXtmStqSn207zxxtSxY4+NVUJjx+auXPmuujpqm/vumzr2xS+qGV9oFJ6Slfr6CMjNm1PHFi+OKUilNCiUrYqKmNK0YEHq2IoV8X/a0JC7csngKTwlI21tsZHH+PGpY0kz/YgjimPC+0gpK4uujfRmfEdHdIGoGZ//FJ4yaM3NsOOOXWtLp56qZvq2UjO+MCk8ZVB6G03fsqX05m4Ol76a8RqNz18KT+lXRwf86lfRlEx2HaioiDd0+gi7bLv0ZnyytDMZjV+0SJst5xuFp/SppSWajhdckDp2yCExWKSt1oZP0oyfNi117Kij4iqhWhufPxSe0qvGxngTr1iROrZoETz9dHFu6JFvxoyB1au7Tqq///7YC1Vr4/ODwlN62LIlgjOp5ZSVRVPy4x/XaPpISibVb9mSOvbhh1Hr13Sm3FN4ykfa2uALX4irQiamTIk3bLHvhJTPamqiJZDeVVJTE1cTbW/PXblKncJTgNQ0pJtuSh27/PJYj13Iu7wXi7FjY0rYscemjl11VVzCRNOZcmNUrgsgudfQEJPe00dzN2/uOhFecm/06FjBtXgx/PVfx7HVq6NVsHmzBvFGmmqeJay/aUgKzvxkFn3PDQ2pTVfa2iJANZ1pZCk8S1RrKxx8cNdpSPvuq2lIhSKZzjRpUuqYpjONLIVnCdq6NWqWzz+fOnbjjbHCRdOQCkdlZVzu+MK0ay/cf3/sdNXUlLtylQqFZ4lJapbpb676+hhl196bhWfUKLjuuq7TmTZv1u5MI0FvlxLR3g7f/37vuyGV6iUyiklNTTTj01sONTXw+OPqBx0uCs8S0NwM06fDlVemjmnT4uJTVRWtiPTdmY4+Or5aW3NXrmKl8CxyyYqUd99NHXviCW1aXKyS3ZnSd/dftChaHFrWObQUnkXKPZps48alNtUtL49+sCOP1DLLYlZWFn3Y6f2gTU2pmqkMDYVnEWptTTXXEtOmaZllqeltWef48VrWOVQUnkUmmYa0aFHq2OWXx0oUXQK49PS1rHOHHTSdaVtpeWYRqa/vuTJoyxZtWlzqelvW+f77Eaz19ZptkS3VPItAe3s0xXqbhqTgFEgt6+xtOtOvfqWLzWVD4VngmpqiCXbVValjmoYkfeltOtMFF8Dee2t3pkwpPAtYbxdlW7RI05Ckf71NZ1q5MgYTGxtzV65Co/AsQOm7ISWS3ZC027sMRjKdqb4+9XpJdmfSqqTBUXgWmKYm2Hln7YYkQ2PcuKhtTpmSOnb00XDggWrGD0ThWUCSZvp776WOaTck2VaVlXHFgMsvTx174QVdM34gCs8CkGzqkd5MHzUqwlS7IclQKC+Hyy6LqW1JMz65Zrya8b3T2y7P9bapx777arWQDI+amlhokX7N+KOPhoMOUjO+uxENTzO70czWm9nLace2M7OFZvZa5/dJ/T1HKamvjyZV+qYeCxaomS7DK7lmfHozftkyNeO7G+ma503ACd2OfRd4xN33BB7pvF/S2trg3HO7TnpPmunnnKNmugy/pBmfPhqfNOO//32tjYcRDk93fxzY2O3wycDNnbdvBuaNZJnyTbIqaMGC1LFDDlEzXXIjGY1Pb8ZfeSVMnaq18flQh9nB3dcBdH6f2teJZnaBmS0xsyV1dXUjVsCRkGwhV13d9UX5xBPw9NNqpkvuVFbCmjVdP9A3bkytjS9V+RCeg+buv3T3WnevnZI+Ma3ANTfD7Nldt5CbNCk+8bX3puSDsrLoMmpo6PpBPn58LMwoxcGkfAjP98xsJ4DO7+tzXJ4RlczdfO211LHLL4e6Oq1Nl/xTXR2v2VNPTR1bvLg0B5PyITzvB87uvH02cF8OyzJiWlpie7CamtQcuvLyaAZddlncFslHo0fDnXf2Pif03HNL57rxIz1V6TbgKWC2mb1lZucBPwGON7PXgOM77xe1pLb5xBOpY4ccoitZSmFJ5oTuvnvq2IIFpVMLHenR9tPcfSd3H+3u0939Bnd/392Pdfc9O793H40vGq2tcNJJ8aJL9k80i2aPBoWkEI0ZE11O6YNJra1RCTjvvOKuhZoX6Lqr2tpaX7JkSa6LMWgNDTBxYtf5cfvuC0uX6vIYUhyamuBjH4O3304dq6yMLRMLadMaM1vq7rUDnZcPfZ5FLb1vMwnOpLb50ksKTikevU1pamqKQaaTTiq+a8crPIeJe/Rpdu/b3Hff6Cc64ghNQZLik0xp6r7N3QMPxHuhoSFnRRtyCs9h0NQEu+4aNc6kb7OsTLVNKR1jx8aeDOm10Pb2aIHNmRNzmwudwnMItbVFJ3lVFaxdmzp+3HGqbUrpSWqhTU1w6KGp4y++GOH6gx8U9hp5DRgNAXd48kn4xCe6rrSoqoL167UmXQSiyb799l3fI2PGxIBSPr1HNGA0QrZuhR137LlEbf78mEScTy8KkVwaNy7mf154YepYc3Mc33PPwttoROGZpZYWOOqoCMf1aQtKDz00XgSXXqpVQiLdjRoF110XlY799ksdX7Uq3ksnn1w4o/IKzwy1tsK8efFpuWhRamnlpEnxqfr00xoQEhlIZWUMnjY0pOaAdnTA/ffHoNIXv5j/E+wVnoPU1hZ/0JoauO++1KdjRQXcdBNs2FBYE4FF8kF1dXRvzZ+faqk1N8P118d77bLL8ndQSeE5gLa2uMzvpEnxB02mWIwaFX03H34IZ5+t3d1FslVeHt1cW7fGbk3JjJSmJrjiCthuu9hpLN9CVG/5PrS0wDHHRGj+6lepyb1m8QdubIy+m1GjclpMkaKR7NbU1BTT+xJJzXT77eGUU/KnT1Th2U1TE+y1V4TmY4+lQnP0aDj++Hj8zjvjvogMvYoKWLgw3muHHZaqiW7eDPfeGyE6d27uJ9orPInmwPz5MGsWTJgAK1ZEzRKi3+Vzn4vm+cMPa+cjkZEyZgw89VSE6PHHxyATxJ63zz8fG+3MmgU335xayTeSSjY8OzrgllvietS77hp9Kq+9lpqruc8+8MMfwgcfwO9+p5qmSK5UVETFpb4evvQlmDkzjjc1xXv2nHNiN6djjoFf/3rkgrSkVhi1t8NVV8Ejj8RlLl55pevjVVWwxx7w7LOpTzkRyS/u8R49/3xYubLn9ZP23juuCTZ3Lnzve5nPtx7sCqOiDs+2Nvj61+N7XV18Sr38ctdzJk6MT61DD4Wrr9YAkEghaW2FM8+E11+Pr02buj6+775xEcUNG+DAAwcXpiURns89t4Rly6KavmJFfCK5w5tvxrG77oqJuN3Nng077xxNcwWmSHFIgnT9eli3LjKhu+OPj93OzKL5bxZfs2fHdMMDDoCyssGFZ0HHxgsvwGc/G4M7mzalVvukSz556uriE2fePDjtNM3LFCk2o0fDHXfE7fZ2+OlPY2Bp++2j1fn44zGKv3Bh158zixZoVRU8+ODgf19Bh+ecObHJal81z/XrVbMUKUXl5XDxxan77hGkDz0U2dBXzXPOnMH/joKOFbPox4DoHBYR6Y1ZZMRQ5oQaryIiWVB4iohkQeEpIpIFhaeISBYUniIiWVB4iohkQeEpIpIFhaeISBYUniIiWVB4iohkQeEpIpIFhaeISBYUniIiWcib8DSzb5rZf5vZy2Z2m5npQhgikrfyIjzNbBrwNaDW3fcDyoHP57ZUIiJ9y4vw7DQKGGtmo4Aq4J0cl0dEpE95EZ7u/jbwf4A1wDpgs7s/3P08M7vAzJaY2ZK6urqRLqaIyEfyIjzNbBJwMrAbsDNQbWZndj/P3X/p7rXuXjtlypSRLqaIyEfyIjyB44C/uHudu7cC9wBH5LhMIiJ9yig8zWzCMJVjDXCYmVWZmQHHAsuH6XeJiGyzTGue75jZDWZ28FAWwt2fAe4Cngde6izXL4fyd4iIDKVMw/NnwPHA02b2584BnHFDURB3v8zd93L3/dz9LHdvHornFREZDhmFp7tfDswETiGmEv2CqI1eZ2YHDHnpRETyVMYDRu7e4e73u/tngD2AfwVOApaa2TNmdo6ZjRnqgoqI5JNtHW3fAmwEGgADJgA3AKvM7OPb+NwiInkrq/A0syPN7BbgbWA+8J/AHHffC9gbeAP49yErpYhInhmVyclm9lXgfxEBuRz4DnCLu9cn57j7SjO7DHhkKAsqIpJPMgpPYgnlvcCX3f2xfs57Dbgi61KJiOS5TMNzhru/N9BJnWvV52dXJBGR/JfpVKUBg1NEpBRk2uf5n/083AFsBpYCNyhoRaSYZdpsN2AWsBPwF+A9YAdiN6R1nfc/DXzTzI5291eGsKwiInkj06lK/ww0AXPdfQ93P8Ld9wAO7jw+H9gTqAOuHNKSiojkkUzD80fA5e7+5/SD7r6UCM4fuftbxBr4o4amiCIi+SfT8JwFbOjjsTrgY523Xweqsy2UiEi+yzQ83wTO7+OxCzofB5gMvJ9dkURE8l+mA0ZXAL8xsxeBu4H1wFTg74D9gNM7zzsOeGaoCikikm8yCk93v83MNhD9mxcDo4FWYAnwSXf/U+ep/wi0D2VBRUTySaY1T9x9IbDQzMqI5vkGd+/odk7TEJVPRCQvDbrP08wqzGyjmZ0EH+3rub57cIqIlIJBh6e7twBtxHxOEZGSlulo+73AqcNREBGRQpJpn+cfgGvM7C4iSNcBnn6Cu/e3/l1EpChkGp53d37/H51fCSfWvTtQPgTlEhHJa5mG598MSylERApMpvM8+9s9XkSkZGQ8zxPAzCYDhwHbAw+4+0YzqwRaNHVJREpBRqPtFn4GvAXcD9wIzOx8+D7gkiEtnYhInsp0qtL3gK8Qa9wPJQaJEg8AJw5RuURE8lqmzfbzgSvc/cdm1n1UfRWwx9AUS0Qkv2Va85wGPN3HYy1oD08RKRGZhufbxNZzvZlDXNdIRKToZRqedwKXmtmRacfczGYB3wJuH7KSiYjksUzD83LgVeBx4LXOY3cCL3Xe/8mQlUxEJI9lOkl+q5kdQ+wY/7fEINH7wA+BW929bchLKCKSh7LZDLkd+HXnl4hIScq02T5szGyimd1lZq+a2XIzOzzXZRIR6UumK4wqzOyyzoBrNLP2bl/b0mz/V+Ahd9+LGLlfvg3PJSIyrDJttv8M+DKxr+c9QPNQFMLMxgNHAefAR7vWtwzFc4uIDIdMw/NU4DJ3v3KIy7E7UAcsMLM5wFLg6+7+4RD/HhGRIZFpn+c44KlhKMco4CDgOnc/EPgQ+G73k8zsAjNbYmZL6urqhqEYIiKDk2l4PkA0r4faW8Bb7v5M5/27iDDtwt1/6e617l47ZcqUYSiGiMjgZNpsvxa4xcw6gP8ANnY/wd3fyLQQ7v6uma01s9nuvgI4Fnhl4J+DZcviuzu8+iq8+SbMmAGrV0NZGVx0EZTrwiAiJa+9HX760/huBjNnRkbMnh33zeCAAwb/fJmGZ9Jkvxy4rI9zso2qrwK3mlkF8AbwhYF+4IUX4LOfhZYWaGuDTZsiRNM9+ijssQesXw+TJ8P778O8eXD66fEfJyLFpaMDfvtb+P3v437SSH3hBXi627ZGZjBxIowaBRUV8OCDg/895t3Tpr+Tzc6h29Uyu3P3mwf/67NXW1vrzz23pM+a51NPwfXXQ2tr7z+/774RpgB77w3XXhv/gSJSWFpb4bTTYMOGCMPGRnj22d7PPfRQ+Mxn+q95lpXZUnevHej3ZhSe/T5R7O85wd17NOWHQ21trS9ZsqTPx93h+efhoYfgrbdSNc9Fi2B5LzNIp02DHXaI27/4BRxySPxnikh+cY9w/NKX4vabb8IHH/Q8b6+9YJ99UjXPUaPg6qsHriSZDS48B6xrmdlG4Dh3f77zvhGX3PhGt/7NWuBJ8uTSw2Ywd258pUv6PR55JG7/5S+wZg28/XZ8ARx2WNRex4yBSy6Bs85SE18kl9rb4Uc/gltvjW661at7njN7Nuy0E2y/PRx00PCPdwxY8+wcHDrM3Z/tvF8OtAK1SaB2Hj8UeNLdRyQ8B6p5DlZ6mG7aBC+91LOpv9NO8Wl1xx0RrKqRigy/jg645Rb4wQ+gqSma5ekqK6PLbfRo+OpXh24cY8hqnsWuvBwuvji+AJqb4ZhjYONGeP31CNd16+KxI46Ipn1FBaxYAWPH5qzYIkWrsRF22y2C8N13uz5WURGPzZgRgzsVFbkpIyg8exgzJgabIJoHJ54YzfoVK+LYe+/F95oa2HFH+NjH4OGHc/tHFCl0zc1w+OExNrFuXdQ60+21F+y/P/zmN1HTzAcKz35UVEQwQqpG+vrrUFcXNdKkn3TcOJgwAe67L14AataLDKyjA26+Gb7zHdi8OaYbpttzz+i/fPTRqNTkm8GG5zQz273zdnnasU1p50wfumLln6RG6g7PPANnnhmd1m1t0Ue6YQMceWTMGdt/f1i4ULVRkd40N8fYweuvQ31918cqK2H6dHjxxfzvFhvsgFH3k6yvY4U2YLQt2tqio/ree3v2zYweHX/8d96Bal1TVEqcOzz5JHz609Gn2b2WudNOcOGFMbsl1ysCh3LAaMCVPqVq1Ci47rr42roV/uqv4I034oXS2hpf48ZF/+jf/A3cdVf+9NeIjISWFjjuuJhz/WG3PdJGj47pRc89FzXOQjNgeI7UiqFCN3YsrFrVe220vh7uvx+qqmJVw0svFeaLRWSwtm6NPst3343xgXT5VMvcFpr6PcSS2ui6dfFJOz2tJ7itLQK2uhrGj4fFi3uuxRcpVB0dcOON0dKqro7B1PTg3GOPaLK/8w5cemlhBycoPIdVVRWsXRsd5Mcdl1oW1tERtdGPfzya9eef37MPSKRQtLbCySdHYJ53HjQ0pCoFlZVw6qnRfF+1Kv8HgTKh8BwBFRUx+t7UBPPnx/SLRGMj3HBDvMjmzImgFSkETU0xz3ns2OiWampKPTZ5Mtx0U7S+7ryzOPv6h2xjkJGWD6Pt26KxMTrL33qr6/GysvgEX7dOo/SSnxoaYrONlpaek9n32CP69Au5hjnY0XbVPHMkadI3NcU2WYmkST9uXNRGb7ih5wtUZKS1t8ca88rK6NNsakq9LsvKoluqubn4mub9UXjm2JgxsUFra2uMQNbUpB5rbo7+0MpKOOmkvvcmFRkuzc2xx2V1dexqlN6ttP320Q3V0lKai0LUbM8z7hGmJ58cy0DTlZfHVKeXX9ZUJxlejY2x+cYHH/Rs+UydGvOZi7VbSc32AmWW2iChsTH6kBLt7bGkbezYqKFqqpMMpY6O1OBldXVcsiY9OPfbL+Zvvvde8QZnJhSeeSyZeJ9MdUrfq7ChIaY6VVdrqpNsm9bW6BaqrIzXUnrTPH2qkRZ3dKXwLADJVKeWluhjGj8+9djWrVFbqKqKFR3p00VE+rN1a1x+ZuxYeOCBrn3qpTDVaFspPAtIeXmszNi8OWqeU6emHmttTY10qkkvfUlWAY0bFx+477zT+yqgujo4+2xdfqY/+q8pUNXV0fe0dWv0RfXWpNcovSRaWuCv/zpeE+ed13WTjlGjSnOq0bZSeBa4ysroi+qtSd/SEs2xqqpohjU25q6ckhsffhiviaoqeOKJrh+kEybAggURmqU41WhbKTyLRHqTvrcNSd5/P2qrNTUaYCp2yVrzmppontfX994037QJzjlHTfNs6b+tCHXfkCR9WklDQ2qAafr0aPZLcWhsjBZGVVWsNW9oSD02aVLxbtCRKwrPIpaM0tfXxy7eU6akHmttjS3DqqtjpcgXv6jaaCFqbYV58+JvmMzNTP87Tp0aIbpxo0bNh5rCswSkT7xP1tIn/Vvu8ca6/vpo5k2Y0HPHb8kv7jGbYsKEaJbfd1/8DRPV1akBIE1oHz4KzxKTrKXfujXm8U2enHqsqQm2bIk35OTJMGuW5o3mk61bYZddYLvtYjbFli3RDE9MmRItjPp6DQCNBIVniSori3l8dXWp2mh6s/799+G11yJIp02Lyy6nv1FlZDQ3w0EHxd9g3LjYwnBT2jVrp05N1TLXr9elr0eSrtsuH9VGIbXP6Lp1MULb3h4Tqd95J96848fDKafEpUZG6dUzLFpb4XOfi6lFvV3PvLIywrTQ980sdKp5ShfpI/Xz58POO6euNdPaGjXS66+P0dvp0yNINQl/2yUDP1OmRF/mffd1HfwZMyYCM1kyqRHz3NOWdDKg/q5PD/Em3m67uCTDww+rr22wmpvhiCOi7/Ltt3tOGzODHXcsjitNFpLBbkmn8JSMtLbCGWfAs89GU757rbOqKoJ07Fi45ZboS1UfXOjogF//OjYVbmqCDRt6DsiVlcWg0Kc/Dddco66RXBhseOpPIxkZPRp+97u43dYGX/sa/OEPqXX2jY2pZaCHHx79p+7R/P/jH0uvVtrcDEceGX3H770XfcndTZ4c/cmf+pQCs5Dk1Z/JzMqBJcDb7n5irssj/Rs1Cn7xi7jd3g4//jHcdls0Q5ML261YEd9Xrox+0lmzoiY6cSI89FDxhWkSlskmwitW9NxToKwsPkwmTICLLoqavJZIFp68arab2T8CtcD4gcJTzfb81twc05uam6N5unZtz3PGj4/r4yTN+r33hmuvLZyaV2srnH561MA3boyVPKtWxYdHd3vuGQE5YwY8+GDxfWgUk4JrtpvZdOAzwJXAP+a4OLKNxoyBp56K2+7RR3r++RGMZrB8eYTM44+nfuaxx2Jy9/77p46ZxYj+6afnrnbW3g4/+QmsWZO6rpR7TEhfv77n+ZWV8UEAEZLXXAMHH6y+32KTN+EJXA38b6BmoBOlsJjFwNFLL6WOJTXTiop4/N13o4n7+uvxle6ee2Lwaffdo7mbHkJmcVG88nL4+7/PPGA7OuD22+Evf4la4erVEYzuMSAG8MILqXmwvTn66Di/oSHWmKtmWRryIjzN7ERgvbsvNbNj+jnvAuACgBkzZoxQ6WQ4pNdMIULstttiOlT69mlr1sDSpVEj7U9ZWTzf2WenwtWsa7eF9yTjAAAMp0lEQVSAOyxblgrHlSvhzTdjGtBAvVdlZbHNWyIJy6uugtpa1SpLUV70eZrZj4GzgDagEhgP3OPuZ/b1M+rzLA1JqC5eHLf7qnk+8wz827/FQFTSZ1pREbXAAw6I+8uWwYknxjLTtrZYvbPddvDNb0YY9lXzLC+HL3wB5s5VSJaCgp3n2Vnz/LYGjCQTHR0xhWrPPTOrec6eDQceqFCUlIIbMBLZFmVl8PnP93+OWQRlYu7c4S2TFLe8C093fxR4NMfFEBHpl6bmiohkQeEpIpIFhaeISBYUniIiWVB4iohkQeEpIpIFhaeISBYUniIiWVB4iohkQeEpIpIFhaeISBYUniIiWVB4iohkQeEpIpIFhaeISBYUniIiWVB4iohkQeEpIpIFhaeISBYUniIiWVB4iohkQeEpIpIFhaeISBYUniIiWVB4iohkQeEpIpIFhaeISBYUniIiWVB4iohkQeEpIpIFhaeISBYUniIiWVB4iohkQeEpIpKFvAhPM9vFzP7LzJab2X+b2ddzXSYRkf6MynUBOrUB33L3582sBlhqZgvd/ZVcF0xEpDd5UfN093Xu/nzn7XpgOTAtt6USEelbXoRnOjObCRwIPJPbkoiI9C2vwtPMxgF3A99w9y29PH6BmS0xsyV1dXUjX0ARkU55E55mNpoIzlvd/Z7eznH3X7p7rbvXTpkyZWQLKCKSJi/C08wMuAFY7u7/nOvyiIgMJC/CEzgSOAv4hJkt6/z6dK4LJSLSl7yYquTuTwCW63KIiAxWvtQ8RUQKisJTRCQLCk8RkSzkRZ+nyFBwh2XL4rsZzJkDL74Y3wFeeAH23z++J+cccEB8F8mUwlPyQns7/OQnsHZt6pgZ7Lxz1/PMYNddYc2a+F5WBrNmwWuvxePf+Q60tEBFBfzTP8HFF8Pdd8djf/d3cNVV8K1vpc752c8iSN1Tv8MdVq+GGTPie/LYO+/ATjvBunWpY7vsAt/9LpSXD8//i+QvhacMq7Y2+NrXInA2bYJJk1KPmcHkybBhQwTT009n/vxmMG4cfPghbL89XHNNhGlS85w9O1XzvPvuqHnOmhXht3JllG3Dhq7hmak77oDDD+/7eT74AD7xiQhyhWzxMN+WV00O1dbW+pIlS3JdjJLX0gKf+lTUHLszi5rk668P7rnKymDePEgWj2Va85w9Gw48cPDNcHf485/h1Vezq3nW1cE9va6F693uu0dNddKk+CBxh/r6CP8pU+C222D06ME/nwwPM1vq7rUDnqfwlP60t0dT909/ijd6UtNraIjvy5bBlh67EPR00EEwfnzfNc8dd4QvfAHmzi2cPkh3WLIkmv6TJvVd81y+PAJ6IFOnwvTpqf9fd6ipidtf/zqceWZ8WMjwUnjKoLW0wAknRG0oCcf6+nhs/Xp4++3+f37sWDjkkJ7Hk3CcNw9OO6103/gdHVGrXLAgujG61zyTGvBA9t47/q+7v2V32w1uv1211qGi8JSPuMOzz8KXvtR7zWjVqlRY9mXqVJg2rWfNs6UFHn0UxowZlqKXhKQGe/HFsHFjz5rn88/H/f5MmgQzZ/Y8XlYGixfr75OJwYanBoyKRGsrfP7z8Je/9Hxs69bBNRtnzYo3bnrNs7wcDj0Urr4aRunVMizM4OCDYeHC3h9vbY0m+/r10UWS/gG4YUP0K3/wQXz1Zvvt42/bm0MOgZ//XH/bbKjmWSDa2uArX4HnnusabklN8KWX+n7zJCoro+nXXUVFjFIffHDh9DdK6OiA3/4Wrr02Qra7l1/u/Xi6adNghx1SoZz0s44bB/vsE89dSuGqZnsBSQZlHnkkdcw9XsA1NXF7+fIY3R3ILrtEP2M6s6h9PPhgBKWUjuZmOOaY+N7d229HbXYgSbiOGxevpWSGgFl8nXsunHFG8fRpKzzzQHs7/PjHsHRpao6je9QQzWDixLi9di288cbgnnPGjOjb6l7zbGqCr34VTj+9eF7EMrza2+GnP40P7WQAC1I1z7q6rosW+rPPPtE1kDxPMig2cWKqNXPQQfC97+X/XFeF5xBLRkyfeKLnnMANG1JTbiZPTtUQX3stmk2DNWNGjJwmz5te81Q4ykjrHq691TxffRXee2/wz7nffrDnnqn7U6b0fP8kYTttGuyxR/Tlj+RrvmTC073rmmWI23femfojrVwZt1eu7Hsy9Jo1XSdFu8dk6GRS9MaN8LvfZVfW9DmOvdU8y8vh2GPhoovy/1NZJF1bG3zjG/DKKzBhQt81zzVrogWWjeOPh49/PPV+TO+X775oorc+e7NYQAGRAbNmRRjPmZPKjPQ9DkomPJct67pmGaKGdsklEVgAmzdDdXVq+se2+MQnun5y9lfzLC/XHEcRSLXc7r2362q0/mqeGzak9iXYFklFBSILJkyAqqrY+yDJjAcfjACN80tkqtKcOV3XLEPc3nXXoa15TpsWy+tGugkhUgzKymJQ6YwzBv8zSeAuXhwr0Iaj5plkRrL/QSYKvuYpIjKUBlvzVB1KRCQLCk8RkSwoPEVEsqDwFBHJgsJTRCQLCk8RkSwoPEVEsqDwFBHJgsJTRCQLCk8RkSwU7PJMM6sDVvfy0GRgwwgXZ1sVWplV3uFXaGUutPJC32Xe1d2nDPTDBRuefTGzJYNZl5pPCq3MKu/wK7QyF1p5YdvLrGa7iEgWFJ4iIlkoxvD8Za4LkIVCK7PKO/wKrcyFVl7YxjIXXZ+niMhIKMaap4jIsCvK8DSzH5rZi2a2zMweNrOdc12m/pjZz8zs1c4y/97MJua6TAMxs8+Z2X+bWYeZ5e0oq5mdYGYrzGyVmX031+UZiJndaGbrzSyD667mjpntYmb/ZWbLO18PX891mfpjZpVm9qyZvdBZ3vlZP1cxNtvNbLy7b+m8/TVgH3e/MMfF6pOZfRL4T3dvM7OfArj7RTkuVr/MbG+gA/h34NvunnfXRDGzcmAlcDzwFvAccJq7v5LTgvXDzI4CGoBb3H2/XJdnIGa2E7CTuz9vZjXAUmBevv4fm5kB1e7eYGajgSeAr7v705k+V1HWPJPg7FQN5PUnhLs/7O5tnXefBqbnsjyD4e7L3X1FrssxgEOAVe7+hru3ALcDJ+e4TP1y98eBjbkux2C5+zp3f77zdj2wHJiW21L1zUND593RnV9Z5UNRhieAmV1pZmuBM4BLc12eDJwL/CHXhSgS04C1afffIo/f2IXOzGYCBwLP5LYk/TOzcjNbBqwHFrp7VuUt2PA0sz+Z2cu9fJ0M4O6XuPsuwK3AV3Jb2oHL23nOJUAbUeacG0yZ81wvF6LN71ZIoTKzccDdwDe6tfzyjru3u/sBRAvvEDPLqnukYK/b7u7HDfLU3wL/D7hsGIszoIHKa2ZnAycCx3qedERn8H+cr94Cdkm7Px14J0dlKVqdfYd3A7e6+z25Ls9gufsmM3sUOAHIeICuYGue/TGzPdPungS8mquyDIaZnQBcBJzk7o25Lk8ReQ7Y08x2M7MK4PPA/TkuU1HpHIC5AVju7v+c6/IMxMymJLNZzGwscBxZ5kOxjrbfDcwmRoNXAxe6+9u5LVXfzGwVMAZ4v/PQ0/k8OwDAzE4BrgWmAJuAZe7+t7ktVU9m9mngaqAcuNHdr8xxkfplZrcBxxA7/rwHXObuN+S0UP0ws48Di4CXiPcbwMXu/h+5K1XfzGx/4Gbi9VAG/M7dr8jquYoxPEVEhltRNttFRIabwlNEJAsKTxGRLCg8RUSyoPAUEcmCwlNEJAsKTxGRLCg8RUSyoPCUomBmn+/cULq5c5PbU8zs0c61y8kmuP/SubFJg5m9a2YPmNleOS66FKiC3RhEJGFmx5HaAOZbxJLRfyX2akz2HB0D1AA/AtYB2wFfAp42s73c/d2RLrcUNi3PlIJnZouBScB+7t7ReexQYmPpx9z9mF5+ppwI1PeAS939X0auxFIM1GyXgtYZggcDdyXBCdC5we2b3c79n2b2jJltIvZN/RAYR2wiI5IRhacUuslE8/y9Xh776JiZfRa4g7hMxOnAoUTo1gGVw19MKTbq85RCtwFoBXbo5bEdiC0JIfbyXOXu5yQPdm7iu91wF1CKk2qeUtDcvZ3Y9PhUM/vo9dzZ5zkz7dQqoqme7ixiX0eRjKnmKcXgMuBh4F4z+3ditH0+kD6C/hAwz8z+BXgQmAt8jdjIWSRjqnlKwXP3PxFXSZ0N3AN8B/gGqWlKAL8CrgT+HngA+AzwWWDziBZWioamKknRSibI9zZVSWRbqeYpIpIFhaeISBbUbBcRyYJqniIiWVB4iohkQeEpIpIFhaeISBYUniIiWVB4iohk4f8DhwTX1oc9uhcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1bd2cc09898>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def RHS(x):\n",
    "    return np.cos(x)-(P/(x))*np.sin(x);\n",
    "\n",
    "## do a scan of K...\n",
    "Kguesses = np.linspace(1e-3,4*np.pi, 10000);\n",
    "band_structure = [];\n",
    "for kguess in Kguesses:\n",
    "    val = RHS(kguess);\n",
    "    if(abs(val) <1):\n",
    "        q = np.arccos(val);\n",
    "        E = kguess**2;\n",
    "        band_structure.append([q,E]);\n",
    "band_structure = np.array(band_structure);\n",
    "\n",
    "plt.figure(figsize = (5,5))\n",
    "alpha = 0.1;\n",
    "plt.plot(band_structure[:,0], alpha*band_structure[:,1], '.b', markersize = 1);\n",
    "plt.plot(-band_structure[:,0], alpha*band_structure[:,1], '.b', markersize = 1);\n",
    "\n",
    "# plt.plot(Kguesses, alpha*Kguesses**2, '.c', markersize = 0.1);\n",
    "# plt.plot(-Kguesses, alpha*Kguesses**2, '.c', markersize = 0.1);\n",
    "\n",
    "# plt.axvline(np.pi, linestyle = '--')\n",
    "# plt.axvline(-np.pi, linestyle = '--')\n",
    "plt.xlabel('qa', fontsize = 16);\n",
    "plt.ylabel('Energy', fontsize = 16)\n",
    "plt.xlim((-np.pi, np.pi))\n",
    "plt.savefig('Konig_Penny_bands.png', dpi = 300)\n",
    "plt.show();\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
